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ABSTRACT 

We carry out a new suite of cosmological radiation hydrodynamic simulations that explores the 
relative impacts on rcionization-epoch star formation of galactic outflows and photoionization heating 
from a self-consistently grown extragalactic ultraviolet ionizing background (EUVB). We compare 
the predictions with observational constraints from the cosmic microwave background, the ultraviolet 
continuum luminosity function of galaxies, and the Lyman-a forest. By itself, an EUVB suppresses 
the luminosity function by less than 50% at z = 6 even if it is orders of magnitude stronger than 
observed. This overproduces the observed galaxy abundance by a factor of 3-5, indicating the need 
for an additional feedback process. We confirm that outflows readily suppress both the EUVB and 
the luminosity function into improved agreement with observations. Population I— II star formation 
can reionize the Universe by z — 6 even in the presence of strong feedback from photohcating and 
outflows. The resulting EUVB suppresses star formation in halos with virial temperatures below 10 5 K 
but has a weaker impact in more massive halos. Nonetheless, halos with virial temperatures below 
I0 5 K contribute up to ~ 50% of all ionizing photons owing to the EUVB's inhomogeneity. Overall, 
star formation rate scales with halo mass Mh as M^~ 1A for halos with Mh = 10 8 ' 2_10 2 M Q . This 
is a steeper dependence than is often assumed in reionization models, boosting the expected power 
spectrum of 21 centimeter fluctuations on large scales. The luminosity function rises steeply to at 
least Afigoo = —13 even in models that treat both outflows and an EUVB, indicating that reionization 
was driven by faint galaxies (M1600 > — 15) that have not yet been observed. Outflows and an EUVB 
interfere with each other's feedback effects in two ways: Outflows weaken the EUVB, limiting Jeans 
suppression of low-mass halos; this leads to overall de-amplification of suppression at early times 
(z > 8). Meanwhile, they amplify each other's impact on more massive halos, leading to overall 
amplification of suppression at later times. Our models cannot simultaneously explain observations 
of galaxies, the cosmic microwave background, and the intergalactic medium. Correcting for dynamic 
range limitations and adjusting our physical treatments will alleviate discrepancies, but observations 
may still require additional physical scalings such as a mass-dependent ionizing escape fraction. 
Subject headings: radiative transfer — galaxies: evolution — galaxies: high-redshift — galaxies: 
photometry — galaxies: star formation — dark ages, reionization, first stars 



1. INTRODUCTION 

Understanding the processes that governed galaxy 
growth when the Universe was less than one billion 
years old represents a central challenge for the upcom- 
ing decade. On the largest scales, the feedback effect of 
galaxies on the temperature, ionization state, and enrich- 
ment of the intergalactic medium (IGM) will be probed 
throu gh absorptio n by neutral hydrogen (iFan et al.l 
I200H and metals (lOhl 120021; lOppenheimer et all 120091) : 
through its imprint on the cosmic microwave back- 
ground (jKomatsu et al.l 1201 lh ; and through redshifte d 
emission by neutral hydrogen (jFurlanetto et al.1 2006). 
On much smaller scales, observations of the Milky Way 
and its environment will decode the signature left b y 
reionization in low-mass objects ([Bullock et al.l [2000h . 
Between these two regimes, direct observations will con- 
strain the formation of the first quasars, gamma-ray 
bursts, and galaxies. 
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With the realization that quasars c ould not have reion- 
ized the Universe by themselves (iMadau et al.1 119991: 
Diiks tra et aT1l2004bl: iWillott et all 120101 iTreister et al.1 
201 ID . much attention is now focused on understand- 
ing the galaxies. Using sensitive measurements from the 
Wide Field Camera 3 aboard the Hubble Space Telescope 
(HST), the Infrared Array Camera (IRAC) aboard the 
Spitzer Space Telescope, and other facilities, a number 
of groups have begun to constrain the abundance and 
colors of continuum-selected galaxies out to z ~ 10 (for 
example, i Finkclstei n et ail 120101: iBouwens et al.l 1201 lbt 
Dunlop et al.l 1201 It [Gonzalez et al.l 120111 : iGrazian et al.l 
20111: iMcLure et all 120111 : lOesch et al.l 120111 ). In a 
complementary approach, several groups have identi- 
fied tens to hundreds of galaxies at z > 6 through 
their bright Lyman-a emission lines. These catalogs 
have already begun to constrain the luminosity function, 
star formation, and clustering properties of faint galax- 
ies (IHu et al.ll201fi lOuchi et al.l l2010l iTilvi et anBoiOl : 
iKashikawa et al.l l2011h . The James Webb Space Tele- 
scope (JWST) will soon push the observational frontier 
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back to even earlier times (G ardner et al.ll2009f ). 

Interpreting these measurements as constraints on 
the feedback processes that regulated early star for- 
mation requires insight from theoretical models. To 
date, models have demonstrated that star formation in 
low-mass halos is influenced by at least two feedback 
processes, namely galactic outflows from star-forming 
regions and photoionization heating (hereafter, pho- 
toheating) owing to an extragalactic ultraviolet back- 
ground (EUVB). Str ong galactic outflo ws are expected 
in low-mass systems (|Dekel fc SiHdl 19861 ) but have in fact 
been obser ved in star-forming galaxies up to the high- 
est masses (|Weiner et al.ll2009h . Cosmological hydrody- 
namic simulations have shown that they are required 
in order to reconcile simulations wit h post-reion i zation 
observations of gala xy abundances (IDave et al.1 120061) 
and scaling relations (IDave et al.ll2006l : iFinlator fe Pavel 
120081: IDave et a l. 2011a ffl) as well as the en richment of 
the IGM (|Qppenheimer fe Da"vl[2006l l2008h . 

Photoheating impacts star format ion by depleting the 
baryon fractions in low- mass halos (iShapiro et al.lll994l; 
Thoul fc Weinberg [1991 IGnedinl 120001: IDiikstra et al l 
2004al lOkamoto et all 120081 ). in fact, the impact of 
an EUVB is often us ed to divide h alos into three cat- 
egories (for example. iHaimanl [20091: Halos with virial 
temperatures in the range 300 < T vir < 10 4 K can use 
molecular hydrogen to cool their gas to star-formation 
densities, but they are susceptible to feedback from a 
background in the H2-dissociating Lyman- Werner bands. 
Halos with 10 4 < T v - lr < 10 5 K can cool their gas through 
collisional excitation of neutral hydrogen, but only in the 
absence of a Lyman continuum background. More mas- 
sive halos can grow galaxies even in the presence of an 
EUVB. For the purposes of this paper, we shall refer to 
halos with T V j r < 10 4 K (virial velocities of < 20 kms -1 ) 
as "minihalos"; halos with T vir = 10 4_5 K (20-64 kms -1 ) 
as "photosensitive" , and more massive halos as "photore- 
sistant" . 

High-resolution radiation hydrodynamic simulations 
now indicate that a single s upernova can unbind a mini- 
halo's entire gas reservoir (iKitavama fc Yoshidal [20051 : 
IGreif et all [20071: IWise fc Abell I2008bl ). The contribu- 
tion of minihalos to cosmological reionization is there- 
fore probably weak, although their contribution to the 
metal enrichment of star- forming gas may be impor- 
tant (Wise & Abel] l2008a| ). Photoresistant halos are 
simpler to model because they respond weakly to an 
EUVB. This means that the largest uncertainty in our 
understanding of reionization involves the amount of 
star formation in the abundant but fragile photosensi- 
tive halos. If the inhomogeneous EUVB does not reach 
a significant fraction of photosensitive halos until late 
times, then the y could have dominate d the star for- 
mation density (Barkana fc Loebl I2000D and hence the 
ionizing photon emis sivity ([Choudhurv fc Ferraral 120071 : 
iMunoz fc Locb 2011) throughout much of the reioniza- 
tion epoch. 

The efficiency of feedback in photosensitive halos is 
also important because of its indirect impact on the more 
massive halos into which they merge. In particular, feed- 
back influences not only the lowest-mass system that is 
permitted to form stars, but also the way in which the 
star formation rate (SFR) scales with the halo mass (Mh) 



at higher masses. This hierarc hical filtering is not ac- 
counted for in idealiz ed models (|Thoul fc Weinberd ll996l : 
Diikst ra~et al.ll2004al) . but it is important because the re- 
sulting SFR-M/j scaling, when modulated by the depen- 
dence of the ionizing escape fraction / esc , determines the 
expected power spectrum of 21-centimeter fluct uations 
at a given ionization state ([McQuinn et al.ll2007[) . 

As such, it is a key ingredient in numerical simulations 
of reionization. Over the past decade, a number of groups 
have used simple models for star formation in dark mat- 
ter halos to model t he reionization of representative cos- 
mological volumes ( Sokasia n et al.l 120031: iGnedin fc Fan! 
2006 ] iKohTer et al.l boOTt iTrac fc Cenl 120071: iTrac et al 



2008T IlIieTet all 120071: 1 Aubertfe Tevssierll2010l: see also 
the review in ITrac fc Gnedinll2009D . These groups have 
achieved impressive success in relating the growth of 
structure across a wide dynamic range to observational 
constraints on reionization from the cosmic microwave 
background (CMB) and the Lyman-a forest. However, 
their predictions remain dependent on the assu med con- 
version from M h to SFR (jMcQuinn et al.ll2007D . and few 
of them have confronted post-reionization observations 
that constrain this scaling such as t he IGM tempera- 
ture or the galaxy LF (see, h owever, iZheng et al.ll2010l 
and iPetkova fc Springelll2010l ) . 

Ideally, models that use galaxies to reionize the Uni- 
verse should be tested against observations of galax- 
ies. For the foreseeable future, the two principal con- 
straints on star-forming galaxies will be the Lyman-a lu- 
minosity function (LF) of Lyman-a emitters (LAEs) and 
the ultraviolet (1350-1600 A) continuum LF of Lyman- 
break galaxies (LBGs). Efforts to date have focused 
more on modeling the LAE LF because its observed evo- 
lution s hould be sensitive to the late stages of reion- 
ization dMalhotra fc Rhoadsl 120061 : iMcQuinn et al.ll2007l : 
llliev et al.l 120081 ). Unfortunately, this exercise depends 
on assumptions regarding both the star formation effi- 
ciency within halos of differe nt masses and the escape 
fraction of Lyman-a photons ([Zheng et al.ll2010l ). Esti- 
mating the contribution of LAEs to reionization depends 
additionally on the unknown escape fraction of ionizing 
photons. In short, it is easier to model the impact of 
incomplete reionization on an assumed intrinsic LAE LF 
than to constrain star formation based on the observed 
one. 

Using the ultraviolet continuum LF (hereafter, the LF) 
to constrain the galaxy contribution to reionization has 
to date been prevented by two concerns. First, the 
LF is an uncertain tracer of star formation owing to 
the unknown amount of dust extinction. Recent work 
has demonstrated that galaxies at z = 6 have quite 
blue UV continua, suggesting that they are relatively 
dust-free dBouwens et aLll2010al : iFinkelstein etafl 120101 : 
iDunlop et al.l 120111 ). If true, then inferring SFRs of 
LBGs from their 1350-1600 A luminosities is less un- 
certain at z > 6 than at lower redshifts. The un- 
certainty in the ionizing emissivity of galaxies then re- 
duces to a single parameter, the ionizing escape frac- 
tion. This parameter remains uncertain, but observa- 
tions now suggest that it wa s larger at early times than 
at prese nt dSiana etajj20nj, and may have been as high 
as 50% ([Rauch et al.ll201ir) . The second obstacle was the 
small available sample sizes at z > 6. Recently, however, 
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observation s have begun to constrain the reionization - 
epoch LF (jMcLure et all 120111: iBouwens et "all 1201 lal ). 
There is now a broad conse nsus as to the abun dance 
of bright galaxies at z — 7 (jMcLure et al.l [201 lh : con- 
straints at z = 8 have emerged; and star-forming 
gal axy candidates have even been identified at z = 
10 (IStark et al.ll2007HBouwens et aT1l2011at lOesch et al.l 
1 2 llf ) - These observations are currently being reinforced 
by even more complete samples from the Cosmic As- 
sembly Near - Infrared Deep Extragalactic Legacy Surve y 
(CANDELS: lGrogin et al.ll2011l . lKoekemoer et al.ll2011f ). 
They do not yet directly constrain star formation in pho- 
tosensitive halos, whose galaxies remain fainter than cur- 
rent detection limits. However, the blue observed UV 
continua indicate that they constrain the star formation 
within photoresistant halos fairly directly, which may be 
invoked as an indirect constraint at fainter luminosities. 

Progress in understanding the significance of photosen- 
sitive halos requires a model for galaxy evolution that 
can be tested against the extensive observational con- 
straints from the post-reionization Universe, and that 
can be extrapolated into the reionization epoch through 
the incorporation of a self-consistently grown inhomo- 
geneous EUVB. Such a model would treat not only 
the spatial inhomogeneity of photoheating, but also 
the no nlinear couplings between different feedback pro- 
cesses (jPieri fe MarteUl2007t IPawlik fc Schavd[2009l) . 

Cosmological hydrodynamic simulations represent a 
mature theoretical model for galaxy evolution that 
can serve as a starting point. Over the last decade, 
numerous studies have shown that the adoption of 
theoretically- and empirically-motivated models for 
galactic outflows brings their predictions into reason- 
able agreeme nt with a wide variety of o bserv ations 
of the IGM (iQppenheimer fc Pavel 120061 I2008D and 
of galaxies dDave et al.l 120061: IFinlator fe Pavel 120081 : 
IFinlator et al.l 120111: iDave et al.ll2011al lbD. These stud- 
ies have generally assume d a spatially-uniform, opt ically 
thin EUVB such as that of lHaardt fe Madaul (|2001l here- 
after HM01). In order to adapt this framework for 
the reionization epoch, we recently developed a time- 
dependent continuum radiative transfer technique that 
is opt imized for cosmological volumes (Fi nlator et al.l 
l2009al) . We used this method to model the growth of 
ionized regions on snapshots extracted from existing sim- 
ulations and verified that outflows leave e nough star for- 
mation to reionize the IGM by z — 6 (Finlat or et al.l 
20093 . We have since integrated this method into our 
custom version of Gadget-2, enabling us to extend our 
previous work into the reionization epoch with improved 
realism. 

In this study, we use this machinery to take a step to- 
ward the assembly of a complete understanding of how 
star formation may have powered reionization by model- 
ing the relative roles of photoheating and outflows within 
the context of three dimensional radiation hydrodynam- 
ics simulations that successfully complete reionization by 
z = 6. We will consider simulations with and without 
each of these feedback processes, for a total of four kinds 
of simulations. This will allow us to explore how strong 
each process is separately as well as how they interact 
when treated simultaneously. 

In § [2 we discuss our custom version of Gadget- 2, 
our radiation transport solver and our approach to mea- 



suring galaxy and halo properties from simulation snap- 
shots. In § [3l we explore how an EUVB and galactic 
outflows impact the baryon mass fractions and star for- 
mation rates of halos, both separately and when consid- 
ered together. We then integrate over all halos in order 
to show how they impact the volume-averaged ionizing 
cmissivity and star formation rate density (SFRD). In 
§ [31 we map these predictions into observable space by 
exploring the normalization and shape of the predicted 
galaxy LF. In § [5l we discuss our simulated reionization 
histories and compare with inferences from the cosmic 
microwave background and the Lyman- a forest. Finally, 
we summarize our results in § [5J We use a standard test 
case to evaluate our code in an appendix. 

2. COSMOLOGICAL RADIATIVE HYDRODYNAMIC 
SIMULATIONS WITH GADGET-2 

In this Section, we describe our numerical methods. 
Briefly, we run a suite of cosmological radiation hy- 
drodynamic simulations that use our custom version 
of Gadget-2 to treat gravity and gas dynamics and 
are coupled with our custom radiation transport solver. 
We identify galaxies in post-processing using skid and 
dark matter halos using fof. Finally, we compute 
the galaxies' broadband photometric colors by convolv- 
ing their star formatio n histo ries and metallicities with 
the lBruzual fc Cha riot (2003) stellar population synthe- 
sis models. We use a tophat filter of width 100 A centered 
at 1600 A to define the UV continuum luminosity. 

2.1. Hydrodynamics, Star Formation and Outflows 

We model structure formation using our custom ver- 
sion of the parallel cosmolo gical galaxy formation code 
Gadget-2 ()Spring cl 2005). This code implements a 
formulation of smoothed particle hydrodynamics (SPH) 
that simultaneously conserves entropy and energy and 
solves for the gravitational potential with a tree-particle- 
mesh algorithm. 

Gas particles undergo radi ative cooling using the pro- 
cesses and rates in Table 1 of lKatz et al.l (|1996f ). We ac- 
count for metal-line coolin g using the collisional ioniza- 
tion equilibrium tables of iSutherland fc Dopital (|1993f) . 
The cooling rate depends on the ionization state, hence 
cooling and ionizations must be computed together. We 
achieve this using a nested subcycling approach that we 
describe below. We initialize the IGM temperature and 
neutral hydrogen fraction to the values appropriate for 
each s imulation's initial redshift as computed by REC- 
FAST (W ong et all 120081 : see E I2.4[) . and we assume that 
helium is initially completely neutral. 

As gas particles cool, they acquire a subgrid two- 
phase interstellar medium consisting of hot gas that con- 
denses via a thermal instability into cold star-forming 
clouds, which are in t urn evaporated back into the hot 
phase by supernovae (jMcKee fc Ostrikerl Il977t ). They 
also acquire the ability to stochastically spawn star par- 
ticles via a Monte Carlo algorithm. This treatment re- 
quires only one physical parameter, the star formation 
timesc ale, whic h is tuned to reproduce th e iKennicuttl 
(|1998l) relation (|Springel fc Hernquistil20i03l) . The phys- 
ical density threshold for star formation is 0.13 cm -3 . 
We account for metal enrichment owing to supernovae 
of Types II and la as well as asymptotic giant branch 
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stars; see lOppenheimer fc Pavel (|2008f) for details of the 
implementation. 

We model galactic outflows (hereafter, "outflows") 
using a Monte Carlo algorithm that applies kicks to 
star-forming SPH particles. We tune the probability 
that particles are kicked and the kick velocity such 
that the resulting outflow mass loading factor and wind 
speeds follow the scalings ex pected for momentum-driven 
winds (|Murrav et al.l [2005). In this model, the ratio 
of the rate at which material enters the outflow to the 
SFR ?7w varies with the host halo velocity dispersion 
a as Vjsss = vo/cr where the outflow amplitude Co = 
150 kms -1 is a free parameter and a is estimated during 
the simulation using an on-the-fly group finder. The out- 
flow velocity is proportional to a. This model broadly 
reprod uces observations of metals in the high-redshift 
IGM (jQppenheimer fc Pavel 120081 : lOppenh eimer et al.l 
1201(1 . 

2.2. Photoionization Feedback 
2.2.1. Discretization 

Our radiative transport (RT) solver discretizes the ra- 
diation field on a three-dimensional Cartesian grid and is 
hence Eulerian. By contrast, Gadget-2 discretizes the 
fluid equations using Smoothed Particle Hydrodynamics 
(SPH) and is hence Lagrangian. We now describe how 
we translate between these two discretizations. 

The conversion from the SPH field to the RT grid oc- 
curs when we compute the emissivity field owing to star- 
forming regions and the opacity field owing to partially 
ionized and neutral gas. The volume- weighted mean 
emissivity 77 in an RT cell is given by an integral over 
the cell's volume V: 

r) = i J 7?(r)dr. (1) 

We derive 77(1*) from the SFRs of star-forming gas parti- 
cles. We account for the metallicity dependence of each 
particle's emissivity using the following fitting function: 

log(Q) = 0.639(- log(Z)) 1/3 + 52.62 - 0.182 (2) 

Here, Q is in s _1 (M©yr _1 ) _1 , Z is the metal mass frac- 
tion, and the last term converts to a lChabrierl (|2003| ) ini- 
tial mass function (IMF). This func t ion re produces the 
equilibrium emissivities of Schaercr (2003;, Table 4) to 
15% throughout the range Z G [10 -7 , 0.04]. We do not 
extrapolate this fit to metallicities outside this range. 
Owing to the fact that gas particles are spatially ex- 
tended and generally straddle cell boundaries, Equa- 
tion[T]becomes, for each cell, a sum over the SPH kernel- 
weighted volume means owing to the SPH particles that 
overlap that cell: 

V = J2 f W(ii - r, hi)r)idV (3) 
. Jv 

Here, r\i is the emissivity of particle i, ri is the particle's 
position, hi is its smoothing length, and W is the SPH 

3 We have verified that doing so has negligible impact on our 
results, which suggests that the ionizing emissivity boost from hy- 
pothetical zero-metallicity stars in atomically-cooled halos has neg- 
ligible impact on cosmological scales. 



smoothing kernel. We evaluate equation [3] by fitting a 
Gaussian to each SPH particle's smoothing kernel so that 
its contribution to each of its neighboring cells reduces 
to a sum over incomplete gamma functions. 

Equations [T}{3] yield the total ionizing photon produc- 
tion rate within a cell, but only a fraction / esc of these 
escape into the IGM. We do not have the spatial resolu- 
tion to model / esc self-consistently, hence we treat it as 
a free parameter and set it to a uniform value of 50% for 
all halos. We will show in Figure [T^] that this choice leads 
to reionization completing by z ~ (7, 6) in our fiducial 
simulation volume (without, with) outflows. 

In our previous calculations, we found that a much 
lower fiducial val ue of ,f esc = 0.12 led to reionization com- 
pleting by z = 6 (Finla tor et al.ll2009b[ ). Two differences 
lead to the higher value that we adopt here. First, our 
post-processing simulations did not account for subgrid 
gas clumping, hence they underestimated the recombina- 
tion rate. In our current simulations, recombinations and 
cooling are computed directly on the SPH particles, au- 
tomatically capturing gas clumping that occurs on scales 
smaller than the radiative transfer grid cells. Second, the 
simulations on which we previously post-processed radia- 
tive transport calculations assumed the HM01 EUVB. 
This background assumes a significant contribution from 
quasars and is harder than our self-consistently derived 
galaxies-only background. The resulting density field, 
which was an input into our post-processing radiation 
transport calculations, possessed a higher temperature, 
significantly more Jeans smoothing, and a lower recom- 
bination rate than occurs in our newer, self-consistent 
calculations. Tuning our simulations to achieve reion- 
ization by z = 6 given our softer background therefore 
requires a higher / esc . 

Each RT cell's opacity \ is also given by a volume- 
weighted mean: 

X=^J v X(r)dV. (4) 

This equation is identical to Equation [T] except that the 
emissivity has been replaced with the opacity, hence we 
evaluate it in the same way. The opacity x(r) is the neu- 
tral hydrogen abundance multiplied by the absorption 
cross section cthi = 6.30 x 10~ 18 cm 2 . We neglect the 
opacity owing to star-forming gas particles because their 
absorptions are implicitly taken into account through 
/ csc . Adopting a smaller cross section in order to match 
the mean energy of ionizing photons (see below) would 
yield a more diffusive EUVB without qualitatively chang- 
ing our results. 

The translation from the RT grid back to the SPH field 
occurs when we compute the photoionization and pho- 
toheating rates. This involves, for each SPH particle, 
summing the contributions owing to the radiation fields 
from each of its neighboring cells. For example, the pho- 
toionization rate Ti for SPH particle i is given by the 
integral (over all space) 

r 4 = J W{vi-v,hi)^)dV, (5) 

where T(r) is the gridded photoionization rate. We eval- 
uate Equation [5] by fitting a Gaussian to each particle's 
smoothing kernel, which again reduces the integral to a 
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sum over incomplete gamma functions. 

The photoheating rate owing to the photoionization of 
hydrogen er HI is the hydrogen photoionization rate Thi 
times the latent heat per photoionization £hi 



CHI 



(6) 



where J v is the mean specific intensity and the integrals 
run from the Lyman Limit j^ll to oo. em depends on 
many factors including the intrinsic slope of the ioniz- 
ing continuum and the imp act of spectral hardenin g in 
the ISM and the IGM (e.g.. iTittlev fc Meiksinl[200l . If 
the simulation volume is large compared to the photon 
mean free path, then all photons are absorbed and we 
may drop the frequency-dependent cross sections <r„ from 
equation |6] to obtain a volume-averaged em- The EUVB 
may be approximated as a power law </„ oc v~ a such 
that chi depends only on the slope a. a can vary be- 
tween 5 for Population I stars and 1.8 for quasars, lead- 
ing to values for em in the range 0.25-0.64 Ryd. By ap- 
plying population synthesis models to the young, metal- 
poor galaxies that dominate our EUVB (§ 12. 3|) . we find 
that their Lyman continua are ch aracterized by a =4— 
5 (see also lBarkana fc Loebll200ll ). This corresponds to 
e H i =0.25-0.32 Ryd. We adopt e H i = 0.3 Ryd, or 4.08 eV 
of thermal heating per ionization. This value agrees well 
with the heating rate per ionization in the galaxies-only 
EUVB derived by H M01, but it is low e r than the value of 
0.6 Ryd assumed bv lFurlanetto fc Oil (2009) or th e rang e 
0.47-2.2 Ryd considered bv iPetkova fc Springell (|2010f ). 
We have not varied em because it is chosen for consis- 
tency with the simulation's predicted stellar populations 
(see § 15.31 for a discussion). 



2.2.2. Solving the Radiative Transport Equation 

We solve for the transport of i onizing radiation usin g 
the moment method presented in lFinlator et a l. (2009a). 
This approach divides the problem into two parts, a mo- 
ment equation update and a long characteristics calcula- 
tion that updates the Eddington tensors. We review this 
technique here and describe our parallelization strategy. 

The zeroth and first angle moments of the radiation 
transport equation in cosmological comoving coordinates 
are: 



87 1 - 

-s-= — V c - T + Airri-cxJ 
at a 



dT_ 

dt : 



-V c ■ (cf J) - c X T 



_ J NfiMVl 



(7) 
(8) 
(9) 



Equation [7] relates the angle-averaged photon number 
density J to the divergence of the photon flux J 7 , the 
angle-averaged photon emissivity rj, and the opacity x> 
where \ accounts for attenuation owing to cosmological 
expansion, redshifting, and absorptions. In our simula- 
tions, r\ depends on the local stellar population and \ is 
dominated by the opacity of neutral hydrogen. Equa- 
tion [5] relates the evolving photon flux to the diver- 
gence of the radiation pressure tensor (the quantity in 
the parentheses) and x- We close the moment hierarchy 
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Fig. 1. — Flowcharts illustrating how the conventional flow of 
a single timestep in Gadget-2 (left) compares with our radiative 
hydrodynamic version (right). 

by relating the radiation pressure tensor to J using the 
Eddington tensor f. We use Equation [5] to derive f from 
the angle-dependent photon number density J\f, which 
we in turn obtain f rom a time-independe nt ray-casting 
procedure (see §3 of lFinlator et al.ir2009a|) . 

We discretize equations [7] [8] fully implicitly in time in 
order to improve the code's stability. We parallelize the 
solution to the resulting linear system of difference equa- 
tions by decomposing the domain of radiative transfer 
grid cells over the set of compute nodes. We then itera- 
tively solve the portion of the linear system that corre- 
sponds to each subdomain and then swap information at 
the subdomain boundaries. We obtain subdomain solu- 
tions from a single Ga uss- Seidel itera tion. This "Addi- 
tive Schwarz Iteration" (Sch warzll870| ) is straightforward 
to implement and lacks global synchronization points 
(because processors only need to swap information with 
processors that host neighboring subdomains). 

The Eddington tensor f encodes the angle depen- 
dence of the radiation field. Errors in the assumed f 
le ad to errors in the sha pes of ionized regions (Figure 5 
of lFinlator elaTl l2009a| ) . In cases where the timescale 
over which the emissivity field changes is long compared 
to the simulation's light-crossing time, obtaining f from 
a time-independent calculation does not introduce large 
errors. The timescale for structure to form is roughly a 
Hubble time, hence this ratio is small for simulation vol- 
umes with comoving side length ~ lOft^Mpc volumes 
at z ~ 10 (note that, throughout this work, we quote 
our simulation volumes' sizes in comoving units). We 
compute f as follows: After each timestep, we evaluate 
whether J has changed within a cell by more than 5%. 
If so, then we update its Eddington tensor by casting 
rays to every source that is not obscured by an optical 
depth greater than 6. We mimic periodic boundaries 
using a single layer of periodic replica volumes. We par- 
allelize this step by first gathering the updated opacity 
and emissivity fields onto all of the compute nodes and 
then assigning roughly equal numbers of updates to each 
compute node. Afterwards, we gather the updated Ed- 
dington tensors onto all of the compute nodes. 

2.2.3. Merging the Radiation and Hydrodynamic Solvers 

Figure [T] illustrates how the addition of our radiation 
transport solver modifies the flow of a single timestep 
in Gadget- 2. In the case of a spatially- uniform, opti- 
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cally thin EUVB (left side), each timestep involves first 
updating the EUVB and then allowing the gas to cool 
and form stars. In the case of radiation hydrodynamic 
simulations (right side), we first update the radiation, 
ionization, and temperature fields simultaneously. This 
step involves three layers of nested loops. The outermost 
loop consists of an iteration between one module that 
advances the gas ionization and thermal state and a sec- 
ond module that advances the radiation field. Within 
the ionization/cooling module, an inner loop alternately 
advances the cooling equations explicitly and then the 
ionization equations implicitly using substeps that are 
limited to 0.002-u/m, where u is the particle's internal en- 
ergy. The ionization solver is itself a Newton-Raphson 
iteration, which we substep separately using a constant 
temperature and ionization rate whenever any of the up- 
dated abundances does not converge to 10~ 4 in 20 itera- 
tions. We substep the outermost iteration whenever the 
photon number density in at least one radiation trans- 
port cell does not converge to 10% within 10 iterations. 
Following convergence of the outermost loop, we proceed 
to its next substep and repeat until all fields have been 
advanced through the timestep. After updating the ion- 
ization, temperature, and radiation fields, we compute 
the number of new "star particles" and "wind particles" 
that are spawned. Finally, we update the Eddington ten- 
sor field. 

Introducing additional physics into Gadget-2 re- 
quires us to slow down the calculation in order to achieve 
a consistent solution. We limit each particle's individual 
timestep based on the evolution of its electron abundance 
as follows: 



\ f< \ 0-ln e (dn e /dt) 1 n e /nj{ > 0.05 
< Omnsidne/dt)- 1 n e /n H < 0.05 



(10) 



Here, n e and uh represent the electron and hydrogen 
number densities, respectively. This criterion evolves 
partially-ionized regions cautiously while limiting the 
ability of neutral regions to slow down the calculation. 
Note that our fully-implicit solution to Equations EHS] 
obliges us to evolve the entire simulation volume's ra- 
diation field whenever any particle's ionization state is 
updated. We also limit the global timestep to be no 
larger than Aln(a) = 0.0035, where a is the cosmologi- 

cal e xpansion factor . 

In iFinlator et al.1 (|2009ah . we demonstrated that our 
radiation transport solver reproduces the analytic solu- 
tion for the growth of ionized regions in expanding me- 
dia, indicating that it conserves photons and accounts 
accurately for cosmological terms. We have verified 
that merging our solver into Gadget-2 preserves this 
agreement. In addition, we have tested our merged 
code's ability to follow the growth of multiple HII re- 
gions in a realistic cosmological density field using Test 
4 from the c osmological radiat ive transfer code compar- 
ison project (jlliev et al.ll2006bl ). We show that it yields 
reasonable agreement with results from other codes in 
the Appendix. 

2.3. Identifying Galaxies and Haloes 

We isolate the galaxies that have formed within 
our volume at each redshift using SKIE0. We infer 

4 http: //www-hpec . astro .Washington. edu/tools/skid. html 



name 


L 1 


winds 


RT grid 


Ml . /M 


time 3 


7/H 4 


r6nWnRT 


6 


no 


_ 


1.4 x 10 8 


1.2 




r6nWwRT16 


6 


no 


16 3 


1.4 x 10 8 


2.1 


3.9 


r6wWnRT 


6 


yes 


- 


1.4 x 10 8 


2.3 


- 


r6wWwRT16 


6 


yes 


16 3 


1.4 x 10 8 


22 


5.1 


r6nWwRT32 


6 


no 


32 3 


1.4 x 10 8 


4.3 


3.6 


r6wWwRT32 


6 


yes 


32 3 


1.4 x 10 8 


31 


4.9 


r3wWwRT32 


3 


yes 


32 3 


1.8 x 10 7 


8.8 


5.3 


r3nWnRT 


3 


no 




1.8 x 10 7 


1.5 




r6nWHM01 


6 


no 




1.4 x 10 8 


1.7 





TABLE 1. Our simulations. All runs use 2 x 256 3 par- 
ticles. The first four simulations explore the relative im- 
pact of ouflows and radiative feedback. The next four 
test our sensitivity to resolution effects. The last run as- 
sumes the optically-thin EUVB of HM01. 

1 in comoving h~ 1 Mpc 

2 virial mass of a halo with 100 dark matter and SPH 
particles. 

3 computation time to z = 6 in 1000 CPU hours 

4 ionizing photons emitted per H atom at xhi,v = 0.01 

the age and metallicity of each galaxy's stellar pop- 
ulation from the age and metallicity distributions of 
its star particles. We compute each star particle's 
contribution to the galaxy's integrated stellar contin- 
uu m by interpolating to its m etallicity and age within 
the iBruzual fc Chariot! (12003[) mod e ls, an d we account 
for dust using the ICalzetti et al.1 (|2000| ) model with 
a n ormalization that i s tied to the stellar metallic- 
ity (jFinlator et al.ll2006D . 

We identify dar k matt er halos following the method 
of IQkamoto et all ( 20081) . First, we identify overdense 
regions using FO!0 with a linking length tied to the 
(redshift-dependent) virial overdensity. Next, we itera- 
tively compute the center of mass and remove the most 
distant dark matter particle until only two particles re- 
main. We then grow a sphere around the midpoint 
between these particles until the enclosed density falls 
to the virial overdensity. Finally, we define the halo's 
baryon mass fraction, SFR, and ionizing luminosity us- 
ing the particles that fall within this radius. 

2.4. Simulations 

Table Q] summarizes our simulation suite. The bulk 
of our discussion will refer to a fiducial cubical volume 
6/i~ 1 Mpc to a side, which we simulate both with/without 
momentum-driven outflows and with/ without a self- 
consistently evolved EUVB. We adopt a cosmology in 
which n M = 0.28, fl A = 0.72, n b = 0.046, h = 0.7, 
as = 0.82, and the index of the primordial power spec- 
tr um n = 0.96. We genera te the initial conditions using 
an lEisenstein fc Hul l) 19991) power spectrum at redshifts 
of z = 249 and 319 for simulations subtending cubical 
volumes of side length 6 and 3 h~ 1 Mpc, respectively. 
We smooth gravitational forces using Plummer equiva- 
lent comoving softening lengths of 469 and 234 h~ x par- 
sees, respectively. The radiative transport simulations 
have been run using Eulerian grids of 16 3 and 32 3 radia- 
tive transport cells. We also ran simulations with eight 
times the mass resolution and one eighth the volume of 
our fiducial 6/i _1 Mpc simulations in order to check nu- 
merical convergence. Finally, we have re-run our fiducial 

5 http: //www-hpec . astro .Washington. edu/tools/f of .html 
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Fig. 2. — The conversion from virial tem perature to halo mass as 
a function of redshift, computed following Barkana & Lockj H1999T ) 
with our simulation's cosmology and a mean molecular weight of 
/J,/m p = 0.59. The upper and lower solid horizontal lines indicate 
the masses of halos containing 100 dark matter and gas particles 
in our 6 and 3 h _1 Mpc simulations, respectively; these represent 
the threshold for halos to possess resolved masses and density pro- 
files (Trcnti et al. 2010). Our fiducial (6/i~ 1 Mpc) volume resolves 
the HI cooling limit for all z < 7. 

simulation without outflows and under the assumption 
of the optically-thin HM01 EUVB for comparison with 
previous work. 

The sixth column of Table Q] indicates that our RT 
solver is fairly efficient. The ratio of the number of 
sources to the number of RT grid cells at z = 6 is 
774/4016 in the r6wWwRT16 and 1526/32768 in the 
r6wWwRT32 simulation. Despite this relatively large 
number of sources, self-consistent RT does not increase 
the total computation time by more than a factor of four 
in the absence of outflows (compare the r6nWnRT and 
r6nWwRT32 simulations). This is because our ioniza- 
tion/RT iteration (the top bubble in the right column of 
Figure [T]) solves the moments of the RT equation for a 
fixed Eddington tensor. Computing the angular depen- 
dence of the ionizing background (that is, updating the 
Eddington tensors) is computationally expensive, but we 
do this for only a fraction of RT cells at each timestep. 

In the presence of outflows, self-consistent RT increases 
the computation time by a total factor of 10-15 because 
outflows generate a significant reservoir of dense gas that 
lives within one optical depth of the nearest source. Its 
brief recombination time increases the IGM's volume- 
averaged recombination rate while forcing the simulation 
to evolve on a shorter timestep. Consequently, outflows 
boost both the computation time and the number of ion- 
izing photons absorbed per hydrogen atom at overlap 
(seventh column in Tabic Q}. 

The goal of the present work is to use a three- 
dimensional numerical model to study the impact of 
galactic outflows and photoheating on star formation 
within photosensitive and photoresistant halos. For ref- 
erence, we show how this mass range varies with redshift 
in Figure [2j The black, solid lines indicate the masses of 



halos that contain 100 dark matter and gas particles in 
our fiducial and high-resolution simulations; halos more 
massive than these thresholds p ossess converged m ass, 
virial radius, and density profiles (|Trenti et al.H2010l ). By 
design, our 6/i _1 Mpc simulations resolve halos with virial 
temperatures above 10 4 K at z < 7. The mass resolution 
limit of our 3/i _1 Mpc volume includes all star formation 
in atomically-cooling halos at all redshifts, although its 
small volume limits the number of massive halos that are 
represented. 

3. OUTFLOWS AND PHOTOHEATING I: THEORETICAL 
INSIGHT 

In this Section, we explore how outflows and photo- 
heating impact star formation. We begin by analyzing 
the baryon mass fraction and SFR as a function of halo 
mass near the overlap epoch. We will find that photo- 
heating primarily impacts halos below 3 x 10 9 M Q , above 
which outflows are the dominant feedback process. We 
will show that outflows and photoheating couple nonlin- 
ear ly at all halo masses, although the sign of the effect 
varies with mass. Next, we weight by halo abundance in 
order to determine how different halos contribute to the 
volume-averaged SFRD and ionizing emissivity. We will 
show that photosensitive halos contribute significantly to 
the SFRD and the ionizing emissivity through the end of 
the reionization epoch. Finally, we study how feedback 
impacts the evolving SFRD. 

3.1. The Halos 
3.1.1. Baryon Fractions and Star Formation Rates 

In Figure [3J we show how the mass fraction in baryons 
/ bar varies with halo mass for each of our simulations at 
two representative redshifts as indicated. The solid red 
curves confirm that all halos retain their baryons in the 
absence of feedback. Allowing a self-consistent EUVB 
to form (heavy red dotted) suppresses the baryon frac- 
tion, with strong suppression visible up to an order of 
magnitude above the HI cooling limit once the universe 
is reionized. The slow growth of f hal with mass owes to 
the hierarchical tendency for halos near the HI cooling 
limit to accrete much of their mass in the form of halos 
that are below the limit. Simulations that treat outflows 
but not an EUVB (heavy blue short-dashed) show that 
outflows suppress baryon fractions only in halos where 
star formation is efficient, hence they are negligible at 
the HI cooling limit. This is an important difference 
between outflows and an EUVB: photoheating primar- 
ily impacts systems near the HI cooling limit while out- 
flows preferentially impact more massive systems owing 
to hierarchical merging. The separate effects of outflows 
and photoheating match at roughly 10 9 Mq, above which 
photoheating weakens while outflows remain efficient by 
design. Finally, treating both outflows and an EUVB 
suppresses the baryon fraction at all sampled masses. 
Curiously, adding galactic outflows to a simulation that 
already grows an EUVB (that is, going from the dotted 
red to the dot-dashed blue curves) boosts the baryon frac- 
tion below 1O 9 M0. We will return to this "suppression 
of suppression" below. 

In order to compare our predictions to previous work, 
we show with a magenta long-dashed curve the baryon 
fraction that results without winds and with the opti- 
cally thin HM01 EUVB. Comparing this with the dotted 
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Fig. 3. — The median mass fraction in baryons (/bar) normalized 
to Q b /QM as a function of halo mass. Red and blue curves indicate 
fiducial simulations with/ without outflows and with/ without a self- 
consistent EUVB as indicated. Thick and thin curves indicate that 
the EUVB is discretized using 16 3 and 32 3 RT cells, respectively. 
Magenta long-dashed curves denote simulations without outflows 
and with an optically-thin HM01 EUVB. Left and right vertical 
dashed lines mark T vir = 10 4 K and 10 5 K from Figure [5] Without 
feedback, halos contain their global baryon fraction. An EUVB 
suppresses the baryon fraction in halos below 3 X 10 9 Mq, with all 
baryons removed from halos a few times more massive than the HI 
cooling limit. Including both feedback processes causes the baryon 
fraction trend to flatten at 0.6Qt/^M f°r halos above IO^Mq. 



red curves shows that a spatially-resolved EUVB yields 
higher baryon fractions at low masses at z = 7. Low 
mass halos have only just been exposed to the EUVB 
in this simulation because it completes overlap around 
z = 7 (Figure I12|) . hence their baryon reservoirs have 
not yet responded. This illustrates the sensitivity of the 
predicted baryon fractions to inhomogeneous reioniza- 
tion scenarios as compared to models that assume a ho- 
mogenous EUVB. By z = 6, the resulting baryon frac- 
tions are within 10% of each other even though the am- 
plitude of the HM01 EUVB is an order of magnitude 
weaker (Figure IT2|) . This confirms that baryon mass 
fractions depend weakly on the EUVB amplitude, as 
opposed to its bias and hardness. This can be under- 
stood as follows: A halo's baryon mass fraction depends 
mostly on its virial temperature relative to the temper- 
ature of the surrounding IGM. The IGM temperature 
depends mostly on the EUVB hardness at the moment 
it was reionized because it decouples from the EUVB af- 
terwards. Therefore, boosting the EUVB amplitude in 
an ionized region does not impact the local baryon frac- 
tions be cause it does not impact th e IGM temperature 
(see also Mesinger fc Diikstral[2008h . 

In Figure SI we show how SFR varies with halo mass 
in the same simulations. In all cases, star formation is 
negligible in halos below the HI cooling limit because 
minihalos cannot cool by collisionally exciting neutral 
hydrogen. Above this cutoff, SFRs rise superlinearly 
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Fig. 4. — SFR as a function of halo mass. Simulations and line 
styles are as in Figure [3] Star formation is suppressed in halos 
that approach the HI cooling limit (left vertical dashed line) even 
though these halos accrete their cosmic complement of baryons. At 
higher masses, the relation between SFR and halo mass approaches 
SFR oc M/j. Overall, halos that are baryon-depleted also show 
suppressed star formation rates as expected. 

over roughly a decade in mass before asymptotically ap- 
proaching a high-mass trend that is generically steeper 
than SFR oc Mh (black long-dashed line segment with 
arbitrary normalization). An EUVB raises the cutoff 
mass roughly a factor of 3 above the HI cooling limit by 
suppressing the baryon fractions in photosensitive halos 
(Figure [3]). Star formation in halos more massive than 
3 x 10 9 M Q is weakly affected by the EUVB. Assuming 
an optically-thin HM01 EUVB and neglecting outflows 
results in SFRs that are ~ 10% higher than predictions 
from our self-consistent EUVB, reflecting slightly higher 
baryon fractions in the presence of a weaker EUVB (Fig- 
ure [12]). Overall, relative SFRs given different feedback 
models follow the trends expected from the baryon frac- 
tions. 

We have fit the simulated trend of SFR versus Mh in 
the wind models both with and without an EUVB using 
a fitting function that consists of a power-law at high 
masses and a turnove r at low masse s (the turnover uses 
the functional form of Gncdin 2000): 
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(11) 



Table [5] gives the parameters of the resulting fits. In ob- 
taining these parameters, we allow the normalization a, 
high- mass slope b, turnover mass M c , and turnover slope 
a to vary independently at each redshift except as in- 
dicated in the caption. The normalization a is 10-20% 
lower in the presence of an EUVB, and it decreases by 
40% from z=10 — > 6. The high-mass slope b is always 
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z 


a 1 


b 


Ml 


a 




r6wWnRT (winds, no 


EUVB) 




10 


0.56 


1.29 


8.08 


3.3 


9 


0.50 3 


1.29 


8.15 


3.6 


8 


0.50 


1.35 


8.18 


5.1 


7 


0.36 


1.29 


8.27 


4.0 


6 


0.29 


1.30 


8.30 


5.1 




r6wW 


ivRT16 (winds + 


EUVB) 




10 


0.49 


1.27 


8.10 


4.6 


9 


0.47 


1.30 


8.18 


4.0 


8 


0.42 


1.36 


8.22 


5.1 


7 


0.31 


1.33 


8.36 


3.0 


6 


0.25 


1.40 


8.61 


3.4 



TABLE 2. Fits using Equation[TT]to the simulated rela- 
tionship between SFR and Mh in models with outflows, 
and with and without an EUVB. 

1 in M Q yr" 1 

2 in log lo (A/ c /M ) 

3 Owing to outliers, we were forced to set a — 
0.50M Q yr _1 by hand in this model. All other parame- 
ters were allowed to vary freely. 

steeper than the linear relation SFR oc Mh- We have 
verified that our higher-resolution winds+EUVB model 
(r3wWwRT32; see Tabled]) also yields b = 1.3-1.4, hence 
this scaling is robust to mass resolution effects. This has 
implications for the typical size of ionized regions and 
the expected power spectrum of 21 centimeter fluctua- 
tions. For example, Figure 17 of IMcQuinn etaD ([2007) 
indicates that changing b from 1 to 5/3 increases the ex- 
pected power by ~ lOx at 0.1 h Mpc -1 scales if the 
neutral fraction is 20%. Our simulations are consistent 
with the steeper end of this range, indicating more power 
at large scales. The cutoff mass M c tracks the HI cool- 
ing mass in the absence of an EUVB, but photoheating 
boosts it by an additional factor of 2 towards z = 6. The 
turnover slope a is not well-constrained because the scat- 
ter grows large to low masses, but it lies within the range 
3-5. These fits may be compared with other models for 
galaxy growth during the reionization epoch. 

Resolution limitations impact our predictions in two 
ways. First, the limited spatial resolution of our radia- 
tion transport solver can dilute the EUVB as long as ion- 
ized regions are not large compared to the grid cells, over- 
suppressing unbiased halos whose environments would 
remain neutral for longer at higher resolution. We show 
this in Figure [3] by using heavy and light curves to com- 
pare results from our fiducial 6/i Mpc volume when the 
RT grid cells span 375 and 187.5 /i _1 kpc, respectively. 
At both redshifts, the no- wind simulation is insensitive 
to spatial resolution because it completes reionization at 
z = 7. By contrast, the simulation that includes outflows 
and an EUVB is just completing reionization at z = 6. 
In this case, halos below 1O 9 M0 retain more baryons at 
higher spatial spatial resolution because their (relatively 
unbiased) environments reionize later. The impact on 
the predicted SFRs is as expected in Figure [4j 

The second resolution limitation involves the limited 
mass and spatial resolution of our hydrodynamics solver. 
Our fiducial simulation misses some star formation at 
redshifts z > 7 because the HI cooling mass falls below 
its 100-particle resolution limit, which delays reioniza- 
tion. For similar reasons, gas clumping is not resolved 
at densities above the limit given by the SPH smooth- 
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Fig. 5. — (top) The impact of mass resolution on / bar at z = 8 
(upper red) and z = 5 (lower blue). Dot-dashed and solid curves re- 
fer to our r6wWwRT16 and r3wWwRT32 simulations, respectively 
(see Table [TJ, while the vertical line segments indicate the 100- 
particle mass resolution limits, (bottom) The impact of resolution 
limitations on SFR. The units of the y-axes are as in Figures [3ti4l 
We have smoothed all curves with a tophat of width 0.2-0.3 mag- 
nitudes. Increasing mass resolution by a factor of 8 suppresses / ba r 
and SFR by a factor of less than < 50%. 

ing length, which delays star formation in all halos. We 
evaluate these limitations using the r3wWwRT32 sim- 
ulation, which trades decreased volume for increased 
mass resolution. Note that this comparison gives an 
incomplete accounting of resolution limitations because 
shrinking the cosmological volume delays the growth of 
structure ow ing to the lack of long -wavelength density 
fluctuations ([Barkana & Locb 2004) while increasing the 
mass resolution accelerates the growth of stellar mass 
owing to the predominance of sm all structures at early 
times ( Springc l fc Hernquisdl2003T ) : these effects can par- 
tially cancel. A more complete accounting will require 
increased dynamic range, which is a goal for future work. 

In the top panel of Figure [5l we show how / bar evolves 
during z = 8 — > 5 in our fiducial (r6wWwRT16, dot- 
dashed) and high-resolution (r3wWwRT32, solid) sim- 
ulations. At z — 8 (upper red curves), increasing the 
mass resolution by a factor of eight reduces the baryon 
fractions by 0.1 dex. That this occurs even though the 
IGM is still essentially neutral in both simulations (Fig- 
ure [T2|) indicates that the dominant effect is the abil- 
ity of simulations with higher mass resolution to resolve 
star formation and drive outflows sooner, drawing down 
their baryon reservoirs. By z = 5, when both simula- 
tions have completed reionization and the fiducial sim- 
ulation has been forming stars for longer, the predicted 
baryon fractions at 10 9 M Q agree well. Below 1O 9 M0, the 
high-resolution simulation predicts higher baryon frac- 
tion be cause its small volume leads to a delayed overlap 
epoch (jBarkana fc Loeb1l2004P ); in essence, halos at the 
HI cooling mass have only just seen an EUVB in the 
high-resolution simulation whereas they have seen one 
for many dynamical times in the fiducial volume. 

In the bottom panel of Figure [SJ we show that the off- 



10 



Finlator et al. 




9.5 
log(M h /M ) 

Fig. 6. — Nonlinear amplification of SFR suppression as a func- 
tion of halo mass in our fiducial volume at three different redshifts 
as indicated. Suppression of the EUVB by outflows weakens sup- 
pression of star formation in photosensitive halos. Coupling be- 
tween outflowing material and an EUVB boosts the suppression of 
star formation in photoresistant halos. 

sets in / bar translate into offset SFRs as expected: At 
z = 8, the high- resolution simulation has a slightly sup- 
pressed SFR while at z = 5 its extrapolation to high 
masses roughly agrees with the trend in the fiducial vol- 
ume. Halos less massive than 1O 9 M0 have significantly 
higher SFRs in the high-resolution simulation because 
their baryon reservoirs have not yet responded fully to 
reionization. In both cases, the SFR vanishes near the 
HI cooling mass. 

Overall, resolution limitations in our fiducial simula- 
tion volume lead to errors of order 0.1 dex in / bar and 
SFR, particularly in halos less massive than 10 9 M Q . The 
important point is that they do not change the quali- 
tative evolutionary trends. We conclude that our fidu- 
cial simulation volume contains enough dynamic range 
to capture the relative impacts of different feedback ef- 
fects even though absolute predictions will have to await 
the arrival of simulations that sample a larger dynamic 
range. 

3.1.2. Nonlinear Coupling of Feedback Mechanisms by Halo 

Mass 

Comparing the blue dot-dashed and red dotted curves 
in Figures [3H4] reveals that accounting for both outflows 
and photoheating leads to higher baryon fractions and 
SFRs below ~ 1O 9 M than in EUVB-only models. This 
suppression of suppression has two causes. First, out- 
flows suppress the EUVB (Figure [1"2")) , leading to a cooler 
IGM with a lower Jeans length that permits lower-mass 
halos to retain their baryons. Second, the simulation is 
just completing reionization at z = 6, hence many of 
the lowest-mass halos have only recently been exposed 
to an EUVB and have already cooled their baryon reser- 
voirs to high densities; the y will continue to for m stars 
for many dynamical times (|Diikstra et al.ll2004al ). Both 
effects weaken the impact of an EUVB on photosensitive 
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Fig. 7. — The differential fraction of ionizing emissivity as a func- 
tion of halo mass in each of our models at z = 8 (top) and z = 6 
(bottom). All curves include the effect of low metallicity in low- 
mass halos, and have been normalized to unit area. Line colors 
and styles are indicated in the legend. Heavy and light curves 
distinguish simulations with 16 3 and 32 3 radiation transport cells, 
respectively. Black vertical dashed lines indicate virial tempera- 
tures of 10 4 K (left) and 10 5 K (right). Vertical segments indicate 
the mass below which 50% of ionizing emissivity occurs. Photosen- 
sitive halos contribute significantly in all models, but particularly 
in the presence of outflows. 

halos at the overlap epoch. They are examples of cou- 
pling between feedback processes, and they illustrate the 
role of radiation hydrodynamic simulations in the study 
of reionization. 

We explore feedback coupling further in Figure[S] Here 
we have computed th e "amplification o f supp ression" 
in two steps following Pawlik fc Schavd (|2009f ): First, 
we divide the mean trend between SFR and halo mass 
from runs that include feedback by the trend from the 
no-feedback run to obtain individual "suppression fac- 
tors". Next, we divide the suppression factor from the 
winds+EUVB model into the product of the suppres- 
sion factors from the EUVB-only and winds-only models. 
This yields the factor by which feedback coupling ampli- 
fies the total suppression of star formation over what 
would be expected if outflows and an EUVB interacted 
linearly. Figure [6] indicates that the sign of feedback am- 
plification depends on halo mass: The tendency of out- 
flows to suppress the EUVB weakens feedback on pho- 
tosensitive halos, while coupling between outflows and 
the EUVB amplifies feedback on photoresistant halos. 
We will return to feedback coupling in § 13.21 We have 
verified that the amplification factors in Figure |6] are es- 
sentially unchanged if we double the spatial resolution of 
the radiation transport solver. 

3.2. Averaging over Cosmological Scales 

3.2.1. The Differential Ionizing Emissivity 

Figures [3]-[4] confirm that photoheating suppresses 
galaxy growth in photosensitive halos. However, they 
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can still contribute significantly to the volume-averaged 
SFRD and ionizing emissivity owing to their abundance. 
We show in Figure [7] the fractional contribution by ha- 
los of different masses to the total ionizing emissivity at 
two representative redshifts. This can be thought of as 
the integral of the product of the curves in Figure H with 
the dark matter halo mass function, but in practice we 
simply sum directly over the simulated halos. We weight 
each halo's SFR by the (self-consistently predicted) metal 
mass fraction of its star-forming gas through Equation^ 

In the absence of outflows (red dotted), an EUVB sys- 
tematically suppresses the total fraction of ionizing pho- 
tons contributed by lower-mass halos, raising the mass 
below which 50% of ionizing emissivity occurs (verti- 
cal segments) from 10 9 ' 7 to 1O 99 M0. The resulting no- 
outflow EUVB is orders of magnitude stronger than the 
optically-thin HM01 EUVB (magenta dashed; see also 
the bottom panel of Figure fl2|) . but it suppresses star 
formation in low-mass halos only marginally more effec- 
tively (as previously seen in § 13.1.11 ) Outflows without 
an EUVB boost the relevance of low-mass halos because 
hierarchical merging systematically depletes the baryon 
reservoirs of more massive halos more effectively (Fig- 
ure [3]) . In simulations that include both forms of feed- 
back, low-mass halos are more suppressed than in winds- 
only models owing to photoheating, but less suppressed 
than in EUVB-only models because the EUVB is weaker. 

We use vertical dashed lines to bound the range of halo 
masses corresponding to T v ; r = 10 4 ~ 5 K, or photosensitive 
halos. Comparing the photosensitive mass range to the 
vertical segments reveals that photosensitive halos con- 
tribute significantly in all models at all epochs, and they 
generate more than half of all ionizing photons for all 
z > 6 in models that include outflows. This fraction is 
slightly overestimated because our fiducial volume does 
not sample the halo mass function above 2 x 1O 1O M0. We 
can estimate the contribution of more massive halos by 
weighting the Sheth-Tormen mass functional by the fits 
in Tableland integrating. Correcting for the limitations 
of our finite volume in this way, we find that photosen- 
sitive halos contribute (49,31)% of all star formation at 
z = (7, 6) in our fiducial winds+EUVB model. This sim- 
ple estimate neglects the weak dependence of metallicity 
on halo mass as well as slight inaccuracies at the mas- 
sive e nd of the Sheth-Tormen mass function (jLukic et al.1 
|2007|) . Nevertheless, it confirms that, for realistic inho- 
mogeneous reionization scenarios, the abundance of pho- 
tosensitive halos more than compensates for their sus- 
ceptibility to outflows and photoheating. 

We have created a galaxy analog to Figure [7] in which 
the x-axis is the absolute magnitude Afi6oo- Summing 
over all simulated galaxies, we find that, in our favored 
winds+EUVB model, half of all ionizing photons origi- 
nate in galaxies fainter than Mi6oo ~ —15. This lumi- 
nosity is roughly ten time s fainter than the deep est cur- 
rent constraints at z = 6 (Bouwcns etalll2011bn . It lies 
close to the resolution limit below which star formation 
is artificially suppressed (Figure ITU)) , hence it is prob- 
ably biased bright. These considerations reinforce the 
need for fainter observations to constrain the dominant 
population of ionizing sources at z > 6. 
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Fig. 8. — (top) The volume-averaged total star formation rate 
density (SFRD) as a function of redshift. (bottom) The ratio of 
the curves in the top panel to the curve from the simulation that 
includes neither outflows nor an EUVB. Line colors and styles are 
indicated in the legend. Heavy and light curves distinguish sim- 
ulations with 16 3 and 32 3 radiation transport cells, respectively. 
Both outflows and an EUVB suppress the SFRD. The solid black 
curve illustrates feedback amplification. The combined suppres- 
sion is less than the sum of their separate effects at early times and 
greater at late times. 

3.2.2. The Volume-Averaged Star Formation Rate Density 

In Figure [8j we show how the volume-averaged star 
formation rate density (SFRD) varies with redshift in 
each of our simulations. The solid red curve in the top 
panel illustrates how the predicted SFRD grows in sim- 
ulations without any feedback. Comparing the dotted 
red and dashed blue curves reveals that outflows and 
photoheating suppress the SFRD by comparable factors. 
The dot-dashed blue curve lies significantly below the 
other curves, indicating that both processes are impor- 
tant in radiation hydrodynamic simulations. Note that 
the predicted SFRD suppression is expected to be robust 
to mass resolution at z < 7 because our simulations re- 
solve the HI cooling limit (Figure [2]). At earlier times, 
they underestimate the impact of photoheating because 
they only partially account for star formation in halos 
below 1.4 x 10 8 M Q . Our high- resolution/small- volume 
simulation has a lower SFRD than the 6/i _1 Mpc volume 
following z — 13 because it lacks sources brighter than 
M 1600 = -15 (Figure EUJ). 

We quantify the amount by which feedback suppresses 
each simulation's SFRD in the bottom panel. By it- 
self, the EUVB suppresses the SFRD by a factor that 
grows to « 30% by z — 6. Its impact flattens below 
z = 7 because photoresistant halos come to dominate 
the SFRD. By contrast, the impact of outflows grows 
with time because the timescale on which galaxies of any 
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mass replace baryons that they eject grows as the cosmic 
density declines. Put differently, winds cause the SFR to 
be dominated by the gas inflow rate at earlier redshifts 
than in models that ignore outflows. By z = 6, outflows 
reduce the SFRD by w 60%. Finally, treating an EUVB 
and outflows simultaneously suppresses the SFRD by an 
amount that grows to sa 75% by z — 6. 

In both panels of Figure |8j we use heavy and light 
curves to compare the results from using coarse and fine 
grids in our radiation transport solver. The resulting 
curves appear coincident in the top panel, but the bottom 
panel clearly shows that lower spatial resolution results in 
decreased suppression of the SFRD at early times owing 
to the tendency for coarse grids to dilute the EUVB. The 
effect is only noticeable while the ionized regions are not 
large compared to the grid cells, however, and by z = 8 it 
is much smaller than the systematic differences between 
different models. 

Our simulati ons predict somewh a t mo re suppression 
by z = 6 than IPetkova fc Springell (|2010t ). who find 0- 
10% depending on their parameter choice (their Figures 
10-11). The difference likely owes predominantly to the 
fact that the mass resolution of our simulations is 37 x as 
high as in their fiducial volume, a choice that reflects our 
emphasis on feedback in low-mass halos rather than the 
IGM. Our high resolution allows us to model the impact 
of an EUVB on the abundant low-mass halos that fall 
below their resolution and that dominate star formation 
at early times (at the expense of subtending a smaller 
cosmological volume). An additional contribution to the 
differe nce may owe to the fact that IPetkova fc Springe! 
(2010) resolve self-shielding filaments more effectively 
than in our simulations, which would weaken the sup- 
pression of low-mass halos in their calculations. 

Figure O suggests that the volume-averaged suppres- 
sion factor from treating outflows and an EUVB simul- 
taneously differs from the sum of their separate impacts. 
This could occur because an EUVB injects entropy into 
inflowing ga s, increasing its cross section to entrainment 
by outflows (jKitavama fc Yoshidal 120051 ) . Alternately, it 
could occur because outflows "puff up" overdense gas, 
increasing its cooling time and rendering it more suscep- 
tible to photoheating. It is expected given Figure [H but 
its overall sign and magnitude depend on the predicted 
SFR-Mh scaling. 

In order to compute the volume-averaged impact of 
feedback coupling, we multiply the EUVB-only (red dot- 
ted) and winds-only (blue dashed) suppression factors 
to obtain the predicted total suppression factor assum- 
ing no feedback amplification (not shown). We then di- 
vide this into the true total suppression factor (blue dot- 
dashed) to find the feedback amplification factor (solid 
black, right y-axis). Overall amplification is below unity 
at early times, indicating that the tendency for outflows 
to suppress the EUVB, thereby boosting star formation 
in low-mass halos, wins over any more direct coupling be- 
tween outflows and the EUVB. This was anticipated in 
Figures [3HH where we showed that, if the EUVB is grown 
sclf-consistently, then photosensitive halos possess higher 
baryon fractions and SFRs with outflows than without. 
The effect is stronger if the radiation transport solver 
uses higher spatial resolution because this further con- 
fines the EUVB to the most overdense regions, weaken- 
ing its impact in voids. Once reionization is well under 



way, amplification grows because star formation is in- 
creasingly dominated by photoresistant halos. Amplifi- 
cation exceeds 20% by z = 6 and continues growing into 
the post-reionization epoch. 

Figures STJ and El may be compared with the results 
of iPawlik fc Schayg (|2009f ). who studied feedback cou- 
pling using a suite of simulations with/without outflows 
and with/without an optically-thin EUVB. They found 
that photoheating and outflows positively amplify each 
other's effects at all redshifts. Our study builds on theirs 
in a number of ways. First, we model an inhomoge- 
neous EUVB, which qualitatively accounts for the pos- 
sibility that the overdense regions that feed gas into ha- 
los could s elf-shield against the EUVB until relatively 
late times (|Finlator et al.l l2009bD . In detail, our simu- 
lations do not yet completely resolve self-shielding be- 
cause our RT cells are wider than the length scales of 
the Lyman Limit sys tems that hos t photosensitive halos 
(~ 10 physical kpc; iSchavd [200 1|) . This explains why 
the baryon fractions in photosensitive halos are not con- 
verged in Figure [3l It will be important in future work 
to model self-shielding within resolution-convergent cal- 
culations. Second, we grow the EUVB self-consistently, 
which accounts for the possibility that outflows suppress 
the EUVB by suppressing star formation (Figure IT2"j) . 
Finally, we account for metal-line cooling. Metal ions 
enhance the cooling rate of enriched gas, weakening its 
response to an EUVB. For the same reason, they boost 
the SFRD and hence the amplitude of a self-consistent 
EUVB. Although we do not isolate the impact of metal- 
line cooling individually, we confirm that it qualitatively 
preserves feedback coupling. 

Quantitatively, our Figure [8] suggests an overall posi- 
tive amplifica tion of 20% by z = 6, as compared to 60% 
in Figure 1 of IPawlik fc Schavd (|2009f ). We attribute the 
difference primarily to the tendency for outflows to sup- 
press the EUVB, leading to "de-amplification" of sup- 
pression in photosensitive halos (Figure [5]) ■ Nonetheless, 
we qualitatively confirm their primary result that hy- 
drodynamic and radiative feedback effects couple nonlin- 
early. Our predicted amplification would probably con- 
verge with theirs at lower redshifts although we have not 
evolved our simulations past z = 5. 

4. OUTFLOWS AND PHOTOHEATING II: 
OBSERVATIONAL IMPLICATIONS 

In this section, we study how outflows and photoheat- 
ing impact the UV continuum LF of galaxies during the 
reionization epoch. We will show that simulations with- 
out outflows overproduce the observed LF while simula- 
tions that include outflows are in reasonable agreement. 
In both cases, an EUVB suppresses the LF normalization 
by less than 30%. Meanwhile, the LF is not expected 
to flatten at luminosities brighter than A/1600 = —13, 
roughly a factor of 100 fainter than current observational 
limits at z = 6. 

4.1. Normalization of the Luminosity Function 

In Figure [9j we show how outflows and photoheating 
impact the predicted LF at two representative redshifts. 
Intuitively, models that predict more strongly suppressed 
baryon fractions (see Figure 13]) also predict more strongly 
suppressed LFs. The no-feedback model predicts the 
highest galaxy abundance within the errors. Next is the 
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Fig. 9. — The simulated rest-frame ultraviolet continuum 
LF including dust for simulations with/ without outflows and 
with/ without an EUVB at two different redshifts (see the legend 
for the line styles). We exclude galaxies with fewer than 64 star 
particles. Simulated errors are Poisson. The top axis indicates 
the integration time for a 5<r detection i n the F200W band o f 
JW ST in units log -i n(t/se c). Data are from|Bouwcns et al.l Q2011br i 
and McLurc ct al. (2010). Both outflows and an EUVB suppress 
the normalization of the LF into improved agreement with obser- 
vations at z < 7. Even in the presence of strong feedback, the LF 
continues to rise steeply below current observational limits. 

EUVB-only model, whose LF is suppressed by < 0.2 dex 
at all sampled luminosities. The third curve corresponds 
to the wind-only model, whose LF is suppressed by a 
factor of ~ 2 at all luminosities. The effect is roughly 
as strong on bright sources partly because hierarchical 
growth spreads the impact on smaller objects into larger 
ones, and partly because we assume that outflows are ac- 
tive at all resolved masses. Finally, including both out- 
flows and an EUVB suppresses the galaxy abundance at 
all luminosities by a total factor of 2-3. 

We also compare these predictions t o recent con- 
straints on the reionizatio n-epoch LF (jMcLure et al.l 
120101: iBouwens et al.l [201 lbfl . At z = 7, the simulation 
including only an EUVB overproduces the observed LF 
in the regime where the simulated and observed ranges 
overlap and the simulation's Poisson errors are not large 
(Migoo ~ —18). This confirms previous indications that 
observations req uire strong hydro dynamic feedback even 
at early times (|Dave et al 1 f2006h . Including outflows 
brings the predicted LF into agreement with observa- 
tions, but now the EUVB's imprint is weak compared to 
current observational errors. The observed LF at z = 7 
has sufficiently large uncertainties that the wind-free sim- 
ulations are only inconsistent with it at the l-3cr level at 
any luminosity, but the fact that the offset is systematic 
argues strongly that it is real. At z = 6, the qualita- 
tive impact of different feedback effects is the same, but 
now observations at the faintest luminosities strongly fa- 
vor models that include outflows. As before, the imprint 
of a realistic EUVB at observable luminosities is weak 
compared to uncertainties. 

The predicted LFs in Figure [9] are affected by the res- 



1 °gio( t 5 y sec ) 

~1> 1 S 1 j 

no winds, no EUVB , 




/ 3 h" 1 Mpc 

,_ fi 6 h~' Mpc 



-15 -10 

^1600 

Fig. 10. — The predicted LF at 2 = 6. Upper solid and dotted 
curves are from our fiducial and high-resolution volumes, respec- 
tively, in simulations without any feedback. Lower long-dashed and 
dot-dashed curves are from simulations of the same volumes that 
include both outflows and an EUVB. For each simulation, heavy 
and light curves indicate the predicted LFs brighter and fainter 
than the luminosity corresponding to our 64 star particle resolu- 
tion limit, respectively. All curves are smoothed with a boxcar 1 
magnitude wide. The top axis indicates the integration time for 
a 5<r detection in the F200W band of JWST in units log 10 (t/sec). 
Resolution limitations boost the luminosities in the fiducial simula- 
tion including feedback (lower long-dashed blue) by a factor of 2—3. 
The predicted LF rises until at least Mi600 = —13, a factor of 100 
fainter than current observational limits at 2 = 6 (Mieoo ~ —18). 

olution effects discussed previously. In order to evaluate 
the strength of these effects, we compare in Figure fTOl the 
predicted LFs from our fiducial and high-resolution sim- 
ulations at 2 = 6. Light and heavy curves show the LFs 
obtained from the complete simulated galaxy catalogs 
and when we exclude simulated galaxies with fewer than 
64 star particles, respectively; we have previously argued 
that this corresponds to the minimum m ass for reason- 
ably converged star formation histories ()Finlator et al.l 
12009) . 

Comparing the lower blue long-dashed and magenta 
dot-dashed curves suggests that the high-resolution 
winds+EUVB simulation yields galaxies that are roughly 
1 magnitude fainter at constant number density than 
the low-resolution simulation with the same physical 
treatments. In other words, halos at z — 6 have 2- 
3 times less star formation at eight times higher mass 
resolution. This numerical effect occurs in the presence 
of momentum-driven winds because higher-resolution 
simulations resolve star formation in lower-mass halos. 
Lower-mass halos have stronger outflows, hence the mas- 
sive halos into which they merge grow systematically 
baryon-deficient (Figure [3]). A qualitatively similar ef- 
fect would occur e ven in the absenc e of feedback (as 
suggested by §3.9 of iKeres et al.l 120091 ) . but momentum- 
driven outflows amplify it significantly. To demonstrate 
this, we also show the LFs from simulations that omit 
feedback (upper curves). In this case, the low- and high- 
resolution simulations agree in the range where they over- 
lap. We conclude that outflows remain required in order 
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to close the gap between the simulated and observed LFs 
at z = 6, although increasing our resolution could allow 
us to decrease the wind normalization a® (§ 12. 1[) . 

In summary, both outflows and an EUVB suppress the 
normalization of the predicted LF without significantly 
changing its shape for absolute magnitudes brighter than 
-^1600 < —15. Over this range in luminosity, the ob- 
served LF favors models that include outflows without 
constraining the EUVB. Resolution limitations artifi- 
cially boost luminosities in our fiducial volume when out- 
flows are included because outflows preferentially evacu- 
ate baryons from low-mass systems. Correcting for this 
numerical effect could enable us to decrease <7o some- 
what, but it would not remove the need for outflows. 

4.2. Turnover of the Luminosity Function 

The LF must flatten below some faint luminosity 
owing to the impact of photoheating on photosensi- 
tive halos and eventually to the HI cooling limit (Fig- 
ure^. Observing this feature would be a major step 
toward constraining the ionizing emissivity that galaxies 
contribute to cosmological reionization because current 
"photon counting" reionization calculations are sensitive 
to the unknown abundance of faint sources (|Chary||2008t 
iMufloz fc Lo^[20TTI: [Bouwens et al.ll2011bO . 

Theoretically, the turnover luminosity depends on the 
reionization history and the nature and stren gth of 
star formation feedback (Bark ana fc Loebl [2000V Re- 
gions that reionize sooner or dev elop a hotter IGM have 
a high er turnover luminosity. Kulkarn i fc Choudhurvl 
(|201 If) considered these effects by combining a semi- 
analytic model for reionization with a Jeans mass cal- 
culation of star formation suppression. They found that 
the mass below which star formation stops varies be- 
tween 1O 8_1O M0 at z = 8 depending on the bias, cor- 
responding to an absolute magnitude between -12 and 
-17. While their results are suggestive, their model had 
difficulty reproducing the evolving LF. Moreover, it did 
not treat gas flows or the growth of ionized regions in 
three dimensions, it did not account for the full temper- 
ature history of a halo's environment, and it assumed 
the timescale s over which halos fo rm stars and respond 
to an EUVB (Dii kstra et al.ll2~004al) rather than comput- 
ing them self-consistently. Our numerical simulations are 
designed to overcome these limitations. 

Observationally, the turnover luminosity at z = 
6 is constrained to b e fainter than -Muy = 
-18 (jBouwens et al.ll2010bD . iMufioz fc Loebl (gop) have 
used a semi-analytic model to map this to a minimum 

n 4 + 0-3 

star- forming halo mass of 10 -»- 9 M Q . This is over an 
order of magnitude more massive than the HI cooling 
limit (Figure [2]). Assuming that some star formation 
continues down to the HI cooling limit, the LF should 
continue to rise to lower luminosities. 

Figure [§] shows that the predicted LF does not turn 
over in any model for luminosities brighter than -M1600 = 
— 16. This constitutes a robust prediction that the ob- 
served LF will continue t o climb to at least this lumi- 
nosity at z — 6 (see also Uaacks et al.l 120111 ). In order 
to trace the predicted LF to even fainter limits, we refer 
to the high-resolution simulation including both outflows 
and an EUVB in Figure[10] This simulation predicts that 
the LF continues its steep rise to at least A/1600 = —13. 



While the light curve flattens at fainter luminosities, the 
LF at luminosities fainter than -13 is dominated by galax- 
ies with fewer than 64 star particles, hence the flattening 
likely owes largely to resolution limitations. 

Figures [SHU confirm that an EUVB suppresses photo- 
sensitive halos, but Figure [TU] indicates that even our 
high-resolution simulation does not fully resolve star- 
forming galaxies at the turnover luminosity. Nonethe- 
less, we may bound the luminosity range within which a 
turnover is expected as follows: The turnover must occur 
in the luminosity range corresponding to T vu = 10 4-5 
K. At z = 6, Equation Qj] and Table [2] indicate that 
the extrapolated (that is, suppression-free) SFRs in ha- 
los with T vlr = 10 4 K and T vir = 10 5 K are 8.3 x 10~ 4 
and O.IOMq yr~ 4 , respectively. The predicted relation- 
ship between M1600 and SFR (M q yr _1 ) at z — 6 is 
Mieoo = (-18.86±0.09) + (-2.56±0.06) log(SFR), hence 
the unsuppressed luminosities would be -11 and -16, re- 
spectively. The bright end of this range may be overly 
conservative given that the LF in our high-resolution sim- 
ulation rises past -16 (Figure H0|). We therefore predict 
that the turnover at z = 6 will occur between -11 and -13. 
Future work incorporating significantly larger dynamic 
range will be required in order to resolve the turnover 
numerically. 

We may draw two conclusions from the faint ends in 
Figures l9l ITOl First, hierarchical formation smooths the 
impact of feedback processes on the LF dramatically such 
that neither outflows nor an EUVB creates a character- 
istic feature above M1600 = —13. Instead, hierarchical 
growth causes any feedback effect that suppresses low- 
mass halos to suppress the normalization of the entire 
LF. We have directly verified that this is true at all 
z < 10. Second, we expect the LF to rise with decreasing 
luminosity until at least M1600 — —13 even in the pres- 
ence of strong outflows and a realistic EUVB. Detecting 
these faint objects at z = 6 will require observations that 
probe a factor of 100 deeper than the deepest observa- 
tions to date, requiring 1-10 Msec on JWST (top axes). 

5. REIONIZATION HISTORIES 

The goal of the present work is to explore the relative 
impact of outflows and a self-consistently modeled EUVB 
on star formation in photosensitive halos. In the process, 
we have simulated cosmological reionization. Our simu- 
lations are not designed to resolve Lyman Limit systems 
or to capture the impact of density fluctuations on length 
scales that are large compared to our simulation volume, 
hence neither the ionization nor the recombination rates 
are expected to be numerically convergent. Nevertheless, 
it is of interest to examine how the resulting reionization 
histories compare to observational constraints. In this 
Section, we compare the volume-averaged neutral hy- 
drogen fraction and EUVB amplitude to observational 
constraints from the Lyman-a forest and the cosmic mi- 
crowave background. 

5.1. Integrated Optical Depth to Electron Scattering 

The probability that photons from the cosmic mi- 
crowave background scatter off of free electrons fol- 
lowing recombination is quantified by the integrated 
optical depth to Thomson scattering, r cs , which is 
a key observable that reionization models must con- 
front. Reconciling the observed cosmic star formation 
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Fig. 11. — The integrated optical depth to electron scattering 
Tes as a function of the age of the universe. Dotted red and dot- 
dashed blue curves indicate radiative transfer simulations without 
and with outflows, respectively. Heavy and light curves compare re- 
sults from simulations with radiative transfer grids with cell widths 
of resolution 375 and 187.5 h~ kpc, respectively. All curves include 
the estimated effect of Helium ionization. Regardless of spatial res- 
olution, helium, and the presence of outflows, simula t ions u nder- 
produce the observed optical depth of Komatsu ct al. (2011). 

history with r es is cha l lenging (IRobertson et al.l 120101 : 
IKulkarni fe Choudhurvl [20lH IHaardt fc Madaul |20T1 
and depends on a number of assumptions including the 
evolving abundance of faint galaxies and the IGM recom- 
bination rate. It is of interest to ask whether our sim- 
ulations reproduce r cs because they model these effects 
self-consistently. 

In Figure llll we show how outflows impact the in- 
tegrated optical depth to electron scattering r cs as a 
function of redshift. All curves assume that helium is 
singly ionized with the same ionization fraction as hy- 
drogen down to z = 3, after which it is doubly ionized. 
Upper red and lower blue curves indicate that simula- 
tions without (with) outflows yield r cs » 0.063(0.054), 
which lies below the observatio nally-determined 68% 
confidence range of 0.073-0.103 (IKomatsu et al.ll20ll . 
Doubling the spatial resolution of our radiation trans- 
port solver (light curves) suppresses r os slightly because 
simulations with higher resolution treat the small self- 
shielded regions (such as Lyman-limit systems) that 
domin ate the IGM opacity at late times more accu- 
rately ([Furlanetto fc Ohll2005D , which slightly boosts the 
neutral hydrogen fraction (Figure [T2"j) . The difference 
is small compared to the discrepancy with observations, 
hence spatial resolution limitations do not dominate the 
apparent underabundance of free electrons at early times. 

Figure [TT] may be compared wit h our post-processing 
calculations (|Finlator et al.|[2009bO . where we also found 
values of r es in the range 0.05-0.07 depending on the 
choice of ionizing escape fraction. That work yielded 
higher values for both r es and the EUVB amplitude while 
assuming a lower / esc because it did not account for sub- 
grid recombinations. Our new calculations treat the IGM 
opacity more correctly, which necessitates a higher / esc 
a nd yields a lower r es and E UVB amplitude. 

iPetkova fc Springel (|2010l ) used a different implemen- 



tation of radiation hydrodynamics into Gadget to sim- 
ulate reionization using 2 x 256 3 particles within a 
10/i _1 Mpc volume. The y accounted for galac t ic out - 
flows using the model of iSpringel fc Hernquistl (2003) , 
hence their simulation is similar to ours, and they found 
t cs = 0.049 (ignoring helium ionization). When we re- 
compute the r es from our r6wWwRT32 model omitting 
helium, we also obtain r os = 0.0493. The remarkable 
agreement is probably coincidental given that our simu- 
lation considered a smaller volume with higher mass res- 
olution, lower latent heat, a different outflow treatment, 
and higher / esc . Nevertheless, it re mains intriguing that 
both s imulations underproduce r es . IPetkova fc Springel 
( 2010) attributed their underestimate partly to the miss- 
ing contribution of ionized helium and partly to the fact 
that sm all cosmological volume s delay the onset of reion- 
ization ■: ISarkaTia fc Loebl 120041 ) . We estimate that he- 
lium contributes Ar cs = 0.005, and our limited volume 
suppresses r os by 003 (assuming a redsh ift delay of 
0.035 from Figure 3 of lBarkana fc Loebl [200l. Summing 
these effects, we estimate that accounting for both he- 
lium and our limited volume would boost the r cs in our 
fiducial volume to 0.058-0.067, still shy of the observed 
range. We will discuss other possible explanations for 
our low T es below. 

5.2. The Neutral Hydrogen Fraction and the Ionizing 
Background 

We show in Figure [12] how the volume-averaged neu- 
tral hydrogen fraction and the ionization rate per hy- 
drogen atom vary with redshift and resolution. If out- 
flows are neglected (heavy red dotted) , then the volume- 
averaged neutral hydrogen fraction drops below 0.01 
around z = 7 while the amplitude of the EUVB at 
z = 6 is overproduced by a factor of 100-1000. Includ- 
ing outflows (heavy blue dot-dashed) delays the comple- 
tion of reionization to z = 6 while bringing both the 
neutral hydrogen fraction and the EUVB into improved 
(though imperfect) agreement with observations. Our 
high-resolution, small-volume simulation (top magenta 
dot-dashed) does not complete reionization until z ~ 5 
owing to its high gas clumping and the lack of long- 
wavelengt h fluctuations. This hist ory is at odds with the 
results of Kashika wa et al.l ()2011l . green arrow). These 
authors found that the Lyman-a LF of LAEs evolves 
more strongly during z = 6.5 — > 5.7 than the UV contin- 
uum LF of LAEs. By assuming that reionization was 
complete at z — 5.7, they used the evolving Lyman- 
a equivalent width to derive a volume-averaged neu- 
tral fraction xhi,v at 2 = 6.5 of ihi,v < 0.4. This 
suggests that our small volume reionizes too late, al- 
though the constraint suffers from uncertainties asso- 
ciated with cosmic variance and the unknown intrinsic 
Lyman-a LF of LAEs. Our small volume also con- 
flicts with the commo nly-accepted reionization redshift 
z = 6 ()Fan et al.ll2006|) . Note, however, that constraints 
from the Lyman-a forest general ly assume a spatially- 
homo geneous EUVB (although see lSrbinovsky fc Wvithd 
120101 ) . This assumption cou l d lead to biased inferences at 
early times. iMcGreer et al.l (|2011[ ) sidestepped the prob- 
lem by obtaining an upper limit to iehi.v from the dark 
pixel fraction in the Lyman-a forest (cyan arrows). This 
maximally-conservative measurement yields much higher 
limits that are in fact consistent with our small volume. 
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Fig. 12. — (top) Volume-averaged neutral hydrogen fraction as a 
function of redshift. (bottom) Ionization rate p er hydrogen atom 
as a function of redshift. Re d squares are fromlFan et all (2006), 
the magenta triangle is from Bolton & Hachnclt (2007), black cir- 
cles are from Srbinovskv & Wvithc (2010), and the (green, cyan, 
black) l imits with (3, 4, 5)-poi nted st ars are from Kashikaw a et afl 
120111 ), IMcGreer etaTl l20Tll). and IBolton et al.l (120111 ) respec- 
tively. We have updated the McGrc er et al.l (120111 ) constraint at 
2 = 5.25 to 0.1 (private communication). The simulated curves' 
colors and styles are as in Figure [TT] The solid black curve indi- 
cates the observed galaxies+QSOs background HM01. Simulations 
without outflows readily reionize the Universe by 2 = 7 while dra- 
matically overproducing the observed EUVB. Simulations includ- 
ing outflows complete reionization around 2 = 6 while improving 
the agreement with the inferred strength of the EUVB. Our high- 
resolution/small- volume simulation completes reionization at 2 = 5 
because it lacks bright sources and resolves more small-scale gas 
clumping. 

However, given that many dark pixels likely correspond 
to regions that are largely ionized, it remains probable 
that our small volume requires a higher / esc in order to 
yield a realistic reionization history. 

Our overall reionization histories could be sensitive to 
the spatial resolution of our radiation transport solver. 
To explore this, we use heavy and light curves to show 
reionization histories from simulations that discretize the 
EUVB using 32 3 and 16 3 grid cells, respectively. Dou- 
bling the spatial resolution of our radiation transport 
solver (light curves) suppresses the EUVB amplitude by 
a factor of *~a few at early times without appreciably 
changing the volume-averaged neutral fraction. We ex- 
pect its impact on our results to be modest. 

It is important to note that, by itself, Figure H"2l yields 
only a joint constraint on the outflow amplitude <7o and 
the ionizing escape fraction / osc because increasing one 
can be compensated by decreasing the other. Break- 
ing this degeneracy requires galaxy observations, which 
constrain the outflow amplitude a®. We have previously 
shown that co — 150 kms -1 suppresses th e predicted 
galaxy LF into agreement with observations (|Dave et al.l 



2006). Figure [TU shows that combining this assump- 
tion with / esc = 50% yields a model that simultaneously 
reproduces observations of galaxies at z =6- 8 and the 
comm only-accepted overlap redshift z = 6 (Fa n et al.1 
2006). These considerations illustrate the crucial role of 
galaxy observations in constraining reionization models. 

5.3. Discussion 

Figure Qj] suggests that the ionized volume fraction is 
too low at early times while Figure [T2l indicates that the 
ionization rate is too high at late times. At the same 
time, we have shown that our current wind model yields 
reasonable agreement w ith the observed UV continuum 
LF o f galaxies at z = 6 (jPave et al.l l2006: Fin lator et al.1 
1 2 1 It Figure |9|), which is proportional to the ionizing 
emissivity. As a guide for future work toward the goal of 
devising a self-consistent model for cosmological reioniza- 
tion and feedback effects, it is instructive to consider how 
our parameter choices could impact our results. The free 
physical parameters are the latent heat per photoioniza- 
tion eui, the ionizing escape fraction / osc , and the outflow 
amplitude oq. The free numerical parameters are the 
simulation volume and the resolution of our SPH and ra- 
diation transport solvers. We now discuss each of these 
in turn. 

Increasing the latent heat per photoionization ehi 
heats the IGM and boosts Jeans smoothing. This si- 
multaneously suppresses the IGM recombination rate 
(and through it the neutral hydrogen fraction) while 
only modestly suppressing the predicted galaxy LF (and 
thr ough it the amplitude of the EUVB). For exam- 
ple, iPetkova fc Springe! (|2010D show that increasing eni 
from 6.4 — s> 30 eV accelerates reionization by Az w 
0.5 while roughly doubling the EUVB amplitude. Our 
current choice eni = 0.3 Ryd is motivated directly 
by the metallicity-dependent stellar continua that our 
simulations predict. The good agreement that simi- 
lar simulations have achieved with the observed his- 
tory of high-redshif t star formation (jPave "elaTl 120061 : 
iFinlator et al J [20 111) and the metal abundanc e in galax- 
ies dFinlator fc Davell2008HDave et al.ll2011b[) and in the 
IGM (Qppcnh eimer fc Pavel I2006D argues that our pre- 
dicted stellar continua are realistic. Increasing eni would 
therefore run against insight from existing observations. 
Moreover, it could lead to implausible IGM tempera- 
tures. For example, IPetkova fc Springel ([2010) find that 
assuming a latent heat per photoionization of 30 eV 
(compared to our 4.08 eV) yields photoionized regions 
with temperatures in excess of 10 6 K, which may be diffi- 
cult to reconcile with the tight observed a ssociation be- 
tween cool gas and Lyma n- Break galaxies ([Weiner et al.l 
l2009t ISteidel et al.1 120101 ) . It would also violate agree- 
ment with the observed distribution of CIV line widths, 
whi ch does not show evidence fo r such significant heat- 
ing (|Qppenheimer fc Pavel 12006=) . Nevertheless, a mod- 
est increase in eni could improve agreement without vi- 
olating existing constraints. Such an increase could be 
associated with spectral filtering through interstellar me- 
dia. 

Increasing / esc causes re ionization to occur so oner by 
strengthening the EUVB (jFinlator et al.1 l2009bl ). This 
means that increasing / csc would boost r os into better 
agreement with observations of the cosmic microwave 
background (Figure [TTJ , but at the cost of exacerbat- 
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ing the disagreement with IGM observations at z ~ 6 
(Figure [T2]) . Decreasing / csc in the case of the no- wind 
model would delay reionization, weaken the EUVB, and 
increase the neutral fraction. This would improve its 
agreement with observations in Figure [T2l but at the ex- 
pense of further suppressing r cs . The only way in which 
/ csc could simultaneously improve the discrepancies in 
Figures lTOJT2l would involve a mass dependence such that 
/csc is larger in low- mass systems whose s tar formation 
halts by z = 6 (similar to llliev "eTall 12001 . 

Weakening outflows (by decreasing er ) accelerates 
reionization in two ways. First, it increases the SFR. As 
can be seen by examining our no- wind simulation — which 
corresponds to the extreme case cr = — this boosts the 
ionizing emissivity. A weaker effect comes from the fact 
that weakening outflows at constant / csc decreases the 
amount of absorbing gas in the IGM. This can be seen 
by comparing the number of ionizing photons per hy- 
drogen at overlap, the last column in Table [1] No- wind 
simulations require w 4 photons per hydrogen to reionize 
whereas wind models require « 5. Intuitively, outflows 
relocate absorbing gas from the ISM into the halo so 
that removing outflows at constant / esc is equivalent to 
boosting / osc . Both effects boost the EUVB amplitude 
and r es . Unfortunately, the EUVB is already too strong, 
and weakening outflows also causes si mulations to over - 
produce the observed LF of galaxies ([Dave et~aT1[2006h . 
which is now well-constrained. In fact, our most recent 
simulations suggest that this parameter is more likely too 
low ( too little suppression) than too high (jFinlator et al.1 
120111) . Hence it is not likely that boosting r es by decreas- 
ing do would lead to overall improved agreement. 

Increasing the simulation volume at constant mass res- 
olution causes reionization to begin sooner by more com- 
pletely sam pling the rare high-q p eaks that were the first 
to collapse (Barkan a fc Loeal2004[ ). This boosts r es with- 
out affecting the galaxy LF or the IGM temperature at 
late times. We estimate that this could increase r cs by 
0.003-0.004, which helps but may not be enough to bring 
predictions and observations into agreement. 

Increasing the mass resolution would help in two ways. 
First, it would boost the IGM recombination rate by 
resolving small self-shielded systems, limitin g the mean 
free path to ionizing photons at late times ([Iliev et al.l 
120051: lAubert fc TeyssieU 12010ft . The potential for im- 
provement can be seen by comparing the light blue and 
magenta dot-dashed curves in the bottom panel of Fig- 
ure [12] These two simulations adopt the same physical 
treatments, but the higher-resolution simulation's (ma- 
genta) increased mass and spatial resolution significantly 
decrease the mean free path to ionizing photons while its 
small volume lacks the sources that dominate star forma- 
tion following z — 13 (Figure [H). The net result is that 
the EUVB amplitude is lower even before its SFRD is un- 
derproduced while the completion of overlap is delayed 
until z — 5. 

The second benefit of increased mass resolution would 
result from higher ionizing emissivity at early times. To 
see this, note that Figures fTTI H2l imply a need for extra 
ionizations at early times but not later on. This could 
be achieved through a population of low-mass halos that 
are acti ve at early time s but suppressed by a mature 
EUVB (|Iliev et all 12007ft . Our fiducial simulations ac- 
count completely for star formation in halos above the 



HI cooling limit at z = 7 of 2 x 1O 8 M , but they fail to 
account for the star formation that must have occurred 
in lower-mass halos at earlier times, when the HI cooling 
limit was lower (Figure [2]) . Increasing our mass resolu- 
tion would boost star formation prior to z = 7 while leav- 
ing it unchanged at later times, thus increasing r es with- 
out affecting the EUVB at z = 6. Unfortunately, for the 
reasons mentioned previously, increasing our mass reso- 
lution would simultaneously increase the ionizing emis- 
sivity and the recombination rate, hence the net effect 
on t cs requires numerical modeling. 

We may estimate how much increasing our dynamic 
range would increase r cs by taking the larger of the 
electron abundances predicted by our 3ft _1 Mpc and 
6/i _1 Mpc simulations at each redshift: 

Ar cs = ca T /(maxK lr ,n e , hr) -n e , r) ^ (12) 

Here, n e ,i r and n e ^ represent the electron abundances in 
the large, low-resolution and small, high-resolution sim- 
ulations, respectively. Carrying out the integral, we find 
that Ar os = 7 x 10~ 6 , with all of the excess scatter- 
ing occurring during the redshift range z =15-22 since 
the smaller volume has a lower SFRD at later times 
(Figure [5]). This correction is small co mpared to the 
discrepancy with iKomatsu et al.l ([20111 ). In detail, it 
is underestimated owing to the decreased presence of 
long-wavelength density fluctuations in our smaller com- 
parison volume. While a full accounting requires sim- 
ulations with increased dynamic range, we conclude for 
the present that increased mass resolution would signifi- 
cantly improve the predicted EUVB amplitude and IGM 
ionization fraction without completely correcting the low 

The final way in which we could reconcile our models 
with observations involves invoking an additional pop- 
ulation of ionizing sources that is active at early times 
such as Population III stars or miniquasars (Mad au et al.l 
2004). In order for miniquasars to bring r cs into agree- 
ment with observations, they must provide an extra 
At os = 0.019. This could be achieved through a constant 
ionized fraction of 0.15 from z = 20 — >• 10. Maintain- 
ing this ionized fraction against recombinations would 
require 2.4 x 10~ 19 ionizations cm _3 s _1 (proper units) 
for a clumping factor of 10, which is a characteristic 
value from our simulations at an ionized volume frac- 
tion of 10%. The ionizing luminosity of Eddington- 
limited accretion by a black hole of mass Mbh is 
~ 10 51 (/ rsr /l) (M BW /1000M^) s" 1 (IMadau et al.ll200l 
Diikstra 2006), hence this level of ionization could be 
supported by a population of 1000M Q black holes with 
a comoving space density of ^1-2 Mpc -3 . This corre- 
sponds roughly to the space density of dark matter halos 
more massive than 5 x 10 8 M Q at z = 6, hence such a 
population could be hosted entirely by atomically-cooled 
halos. 

In summary, Figures E] [12] indicate that our new sim- 
ulations underproduce the ionization rate at early times 
while overproducing it at late times. Increasing eni, the 
comoving volume, or the mass resolution would yield 
simultaneous improvement in both respects. However, 
self-consistency argues against significantly higher values 
of em; simple estimates indicate that the impact of our 
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limited volume is small; and a high-resolution simulation 
suggests that increasing the mass resolution boosts the 
IGM recombination rate enough to counter the boosted 
ionizing emissivity at early times. Meanwhile, increasing 
/ csc or decreasing <tq would improve r es while degrading 
agreement with observations of galaxies and the IGM at 
z = 6. We conclude that simultaneously matching obser- 
vations of galaxies, the IGM, and the cosmic microwave 
background may require us to modify assumptions re- 
garding feedback, the nature of the ionizing sources, or 
the ionizing escape fraction. We plan to consider these 
possibilities in future work. The important point for 
our present purposes is that our simulations yield suf- 
ficiently realistic reionization histories that we may use 
them to gain qualitative insight into how outflows and 
photoheating modulated galaxy growth during the reion- 
ization epoch. 

6. SUMMARY 

We have used a suite of cosmological radiation hydro- 
dynamic simulations to study the impact of galactic out- 
flows and photoheating by a self-consistent, spatially- 
inhomogeneous EUVB on star formation during the 
reionization epoch. The major improvements of our work 
ove r previous rad i ation hydrodynamic studi es (for exam- 
ple. lGnedinll2000t iPetkova fc Springelll2010H are that our 
model includes a treatment for galactic outflows that has 
previously been tested extensively against observations of 
galaxies and the IGM from z = — > 7 (§ [T), and that 
it essentially resolves the HI cooling limit at all relevant 
redshifts. It is the first study to demonstrate agreement 
with the observed UV continuum LF, achieve reioniza- 
tion by z = 6, and resolve the majority of the star for- 
mation that occurs in photosensitive halos. Our major 
results are as follows: 

• By z = 6, a self-consistent, inhomogeneous EUVB 
significantly suppresses the baryon reservoirs and 
SFRs in halos less massive than 3 x 10 9 M Q . Mean- 
while, the star formation activity in halos more 
massive than 3 x 10 9 M Q , which host currently ob- 
servable galaxies, is essentially unaffected. This 
suggests that an EUVB acting alone cannot sup- 
press the UV LF into agreement with observations 
at z = 6. 

• Momentum-driven outflows can suppress the 
EUVB and the LF into improved agreement with 
observations without preventing Population I— II 
star formation from reionizing the Universe by 
z = 6 as long as we assume that, on average, 50% of 
ionizing photons escape from the low-mass galaxies 
that drive reionization. 

• Even in the presence of strong outflows, an inho- 
mogeneous EUVB allows photosensitive halos to 
contribute up to ~ 50% of all ionizing photons 
throughout the reionization epoch. 

• For halos in the mass range M/, = 10 8 - 2_10 - 2 Mq, 
SFR scales as M^ Z ~ 1A . If the escape fraction is 
constant, this implies substantially more power in 
21-centimeter fluctuations on large scales (0.1 h 
Mpc" 1 ) than would be expected if SFR oc M^°. 



• High-resolution simulations indicate that the LF 
continues to rise steeply until at least Mi 60 o = —13. 
This implies that reionization was dominated by an 
abundant population of faint galaxies that has not 
yet been observed. 

• Outflows and an EUVB couple nonlinearly in two 
ways: First, outflows weaken the EUVB, boosting 
star formation in photosensitive halos and lead- 
ing to overall de-amplification of suppression at 
early times. Second, they amplify the impact of a 
gi ven EUVB on photoresis tant halos as suggested 
by iPawlik fc Schavd (|2009() , leading to overall am- 
plified suppression of the cosmic comoving SFRD 
near the end of reionization. By z = 6, overall 
suppression of star formation is 20% greater than 
would be expected from the two effects acting sep- 
arately. 

• Simulations that achieve reionization at z =6-7 
(that is, our fiducial volume) generically underpre- 
dict the optical depth to Thomson scattering while 
slightly overproducing the inferred amplitude of the 
EUVB and the volume-averaged ionized hydrogen 
fraction at z = 6. 

• Discrepancies between observations and models 
could be alleviated through adjustments to phys- 
ical or numerical parameters. In general, how- 
ever, our physical parameters are either indirectly 
constrained observationally or could only improve 
agreement significantly with one observable at the 
expense of another. Simple estimates suggest that 
improved dynamic range would not boost the pre- 
dicted r es into agreement with observations al- 
though it would help. 

• Observations thus suggest the need for an addi- 
tional physical scaling that preferentially boosts 
the ionizing emissivity at early times such as a 
mass-dependent ionizing escape fraction. Alterna- 
tively, an additional population of ionizing sources 
such as miniquasars could have provided the extra 
ionizations. 

While our results represent a significant step forward 
in the modeling of reionization-epoch star formation, our 
predictions are not yet completely converged with respect 
to mass, spatial, or spectral resolution. The need for 
higher mass resolution can be seen in our predicted UV 
luminosity function, which is slightly high owing to the 
delayed onset of star formation at finite mass resolution 
(Figure [To]) . We have also argued that improving our 
mass resolution would alleviate tensions with inferences 
from the cosmic microwave background f £|5.3p . 

The need for improved spatial resolution within our 
radiation transport solver can likewise be seen in two 
ways. First, comparing the number of photons re- 
quired to achieve reionization in our r6wWwRT16 and 
r6wWwRT32 simulations (Table [IJ column 6) reveals 
that the higher-resolution simulation "consumes" 4% 
fewer photons. This is because the low-resolution simula- 
tion smooths over regions that are in reality self-shielded, 
thereby overestimating gas clumping. Second, the fact 
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that the light and heavy blue dot-dashed curves in Fig- 
ures [3]-[4] are not coincident in the photosensitive range 
shows that our simulations do not completely resolve the 
self-shielded regions that host photosensitive halos. A 
convergent prediction would require RT cell widths that 
are smaller than the lengt h scales of Ly man Limit sys- 
tems (~ 10 physical kpc; iSchavd [20011 ). which is four 
times smaller than the RT cells in our r6wWwRT32 sim- 
ulation. Current computational resources already allow 
us to augment our dynamic range somewhat. This will 
improve our estimate of the IGM recombination rate, 
confine the EUVB more completely to overdense regions 
at early times, and resolve the earliest stages of star for- 
mation in photosensitive halos. 

Another potentially important effect involves fluctua- 
tions in the ionizing background on still smaller scales. 
For example, the interstellar radiation field within an 
individual galaxy can dominate the EU VB in wave- 
lengths to which its ISM is optically thin (Gnedinl l2010l : 
ICantalupol 12010ft , but our current simulations do not 
treat it separately from the EUVB. Consequently, they 
may overestimate the amount of feedback that is required 
in the form of outflows in order to reproduce the ob- 
served LF. This process would also assist in alleviating 
lingering tensions between the predicted and observed 
gas fractions of low-mass galaxies, which re main too low 
at low redshifts in models similar to this one (|Dave et al.l 
l2011a| ). We will consider in future work whether this ef- 
fect can be incorporated into our subresolution treatment 
for star formation. 

Upgrades to our radiation transport solver allow- 
ing multifrequency calculations and the addition of 
a background from quasars will result in a more 
accurately- modeled EUVB. Accounting for the longer 
mean free path of higher-energy photons would raise 
the IGM te mperature in ionized as well as neu- 
tral regions (lAbel fc Haehnelt! [19991 llliev et al.l l2006bl: 
iTittlev &: Meiksinll2007D . modifying the relative roles of 
photosensitive and photoresistant halos. 



Finally, recent observations suggest that the way in 
which dense gas collapses into stars could include depen- 
dencies such as the gas metall icity (|Gnedin fc Kravtsovl 
l2010t | Krurnhofz fc D ekcl 20lT[) and the stellar mass den- 
sity (jShi et al.l 1201 lh . Incorporating these dependencies 
could modify the predicted SFR-M^ relation by delay- 
ing star formation activity in photosensitive halos, which 
would again boost the expected power spectrum of 21 
centimeter fluctuations at large scales. 
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APPENDIX 
TEST OF MULTIPLE SOURCES 

Test 4 of llliev et al.1 (|2006b| ) involves evolving the ionization fronts from 16 sources embedded in a static cosmological 
density field, hence it tests the code's ability to resolve the impact of density fluctuations on the shape of the ionization 
front. It also tests the code's treatment for photoheating and radiative cooling. For multifrequency techniques it tests 
spectral filtering, but as our code is currently monochromatic, it illustrates the s ystematics associated wi th neglecting 
the ability of high-energy photons to preheat neutral regions (|Iliev et al ■l l2006bt iTittlev fe Meiksin] l2l)07l) . 

We represent the density field using a uniform grid of SPH particles whose masses vary in such a way as to reproduce 
the test requirement in the absence of SPH smoothing. During run-time, our code extracts the gridded densities from 
the SPH particles using a "clouds-in-cells" approach. This smooths the density field on a length scale of ~ 3 times the 
mean particle separation. Consequently, the opacity field that our radiation transport solver "sees" does not exactly 
reflect the test requirement, and we do not expect our results to agree quantitatively with those of other codes. We 
have verified that increasing the number of SPH particles per grid cell from 1 to 8 yields ionized regions whose shapes 
agree more closely with reference results. We enclose the test volume in 4 layers of opaque boundary cells that prevent 
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Fig. 13. — Ionization state (left) and temperature (right) of the Test 4 case at t = 0.05 Myr (top) and t = 0.2 Myr (bottom). This 
test w as run with 256 3 SPH particles and 1 28 3 RT cells. The ion ization front is significantly smoother than those produced by other 
codes (Ilicv ct al. 2006b; Wise & Abel 2011; Pawlik & Schayc 2011) owing to the diffusive nature of our radiation transport solver, and 
the IGM behind the ionization front remains cold because our monochromatic solver does not treat spectral hardening. Nonetheless, the 
volume and temperature of the ionized region are reasonable, indicating that the code accurately evolves the IGM's thermal and ionization 
evolution. 

photons from "wrapping around" because our code is optimized for periodic boundary conditions. The latent heat per 
photoionization is 1.177 Ryd, which is appropriate for a T = 10 5 K blackbody in the optically-thick approximation. We 
smooth the Eddington tensor with a 27-cell tophat filter in order to suppress rapid spatial fluctuations in the radiation 
pressure tensor. This step ensures photon conservation at the cost of rendering the radiation field more diffusive. We 
have found empirically that smoothing is not required in our current cosmological simulations because their spatial 
resolution is generically much lower than that of Test 4, preventing the appearance of sharp peaks in the radiation 
pressure tensor. However, future cosmological simulations at higher spatial resolution may require smoothing. We run 
the test in an expanding Universe at z = 8.84922 with Qb = &m = 1, and we adjust Hq and the box size so that the 
proper density and volume match the required test conditions. Strictly speaking, Test 4 is meant to be a test of static 
density fields and should not be run in an expanding frame. However, its duration of 0.4 Myr is small compared to 
the Hubble time at this redshift, hence cosmological effects are negligible. 

We show maps of the neutral hydrogen fraction and temperature in a slice of the simulation volume in Figure 1131 
The i onization maps (left column) show similar morphologies to the results from other codes (see, for example, Figure 
19 of Wis e &: Abell l201ll or Figure 10 of [Pawlik fc Schavell201lD . This broad agreement indicates that our approach 
accounts reasonably well for fluctuations in the EUVB on length scales that are large compared to individual cells. 
In detail, our approach yields smoother ionization fronts than ray-tracing codes partly because SPH smooths the 
density field on scales comparable to the mean particle separation, and partly because the moment method smooths 
the radiation field on scales comparable to the cell length. While this motivates further work in order to capture 
fluctuations at length scales that are closer to the size of the RT grid cells, we do not expect it to alter the interaction 
between galaxies and the EUVB when averaged over cosmological volumes. 

For a more quantitative comparison, we show in Figure [TJ] how the volume-averaged and mass-average d ionized 
fracti ons evolve in our own test (black solid) and in tests conducted using enzo+moray ( Wi se fc Abell [20111 ) and C 2 - 
RAY (jMellema et al.l [20061 ). Two kinds of differences are apparent. First, the predicted reionization topologies differ. 
The ratio Xhii,m/^hii,v exceeds unity in a ll codes at early times, indicating that the initial stages of reionization 
are robustly "inside-out" (jlliev et al.ll2006aD . The ray-tracing codes switch to an "outside- in" topology by the time 
the ionized fraction reaches 60%. This indicates that ionization fronts have ar r ived in underdense regions and is 
reminiscent of the "inside-outside- middle" topology discussed in iFinlator et al.l (|2009bD . Meanwhile, our moment 
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Time/kyr 

Fig. 14. — Mass- weighted (light) and volume-weighted (heavy) ionized fraction versus time in our approach (solid black) and using C 2 -RAY 
(dashed red) and ENZO+MORAY (dotted blue). Solid black and dot-dashed magenta curves indicate test runs with 1 and 8 SPH particles 
per grid cell, respectively. The topology of reionization and the ionized fractions are both sensitive to the density field's accuracy. Broadly, 
our predicted ionized fractions are somewhat lower than in the reference codes, indicating that our method predicts a slower reionization 
history in this test case. 
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Fig. 15. — Histograms of neutral hydrogen fraction (left) and temperature (right) at 0.2 Myr. Colors and line styles are as in Figure [141 
Reionization occurs somewhat more slowly in our code as compared to the other codes, and the absence of high-energy photons leads to a 
systematically cooler IGM. 



approach maintains a purely inside-out topology until the ionized fraction reaches nearly 70%, indicating significantly 
stronger ionization-front trapping in overdense regions. Increasing the number of SPH particles per grid cell moves 
the transition to earlier times, suggesting that the difference between the predicted topologies may owe largely to our 
smoothed density field. The second difference involves the reionization rate. Our method predicts that reionization 
proceeds more slowly than in the reference codes, with the neutral fractions ~ 10% higher at the end of the test. 
Comparing the solid black and dot-dashed magenta curves reveals that increasing the accuracy with which we capture 
the required density field weakens ionization-front trapping and accelerates reionization. However, the difference 
is slight, suggesting that increasing our spatial resolution indefinitely would not necessarily bring our results into 
agreement with the reference runs. 

In order to explore our code's slower reionization rate further, we show histograms of neutral fraction and temperature 
at 0.2 Myr in Figure [151 The distribution of neutral fractions has a similar shape in all four cases, but our code shows 
fewer cells in all bins of ihi.v except the neutral bin xhi,v =0.98-1.0. In fact, the fraction of cells in which shi.v = 1 
is ~ 40% in our test runs and < 10 -3 in both reference runs. This dominates the difference between our code and 
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the reference codes in Figure 1141 Accounting more completely for partial ionization by high-energy photons in neutral 
regions would improve agreement. For similar reasons, our code does not pre-heat neutral regions, leading to an 
enhanced volume fraction at low temperature with respect to the reference codes. 

In summary, our radiation hydrodynamic code yields reionization histories that are similar to results from other 
codes. This implies that our approach accounts reasonably well for fluctuations in the EUVB on scales larger than an 
individual cell. In detail, however, underdense regions are colder and more neutral in our code than in others. This 
partly reflects the difficulty of translating a gridded density field into a grid of SPH particles. Mostly, however, it 
simply reflects the absence of high-energy photons that broaden ionization fronts and preheat neutral regions. 



